# R code for paper "Reserves Were Not So Ample After All" by Adam Copeland, Darrell Duffie, Yilin (David) Yang
# Comment time stamp 09/26/2024
# If you have comments or questions, reach out to Yilin (David) Yang at yiliyang@cityu.edu.hk or his latest email address. 

# To replicate the tables and figure, just run every line of this code. You may need to install some packages. 
# Codes for Table 2 Column 2, Figure 7, Figure 8 and Table 14 are included in other stata files. 
 
# Internet Appendix Figure 12, Figure 13 Figure 14 and Figure 15 use confidential data and cannot be directly replicated. 

Fileaddress = dirname(rstudioapi::getActiveDocumentContext()$path)
setwd(Fileaddress)
graphics.off()
rm(list = ls())
gc()
library(Formula)

library(clipr)
library(expss)
library(laeken)
library(rlang)
library(ggplot2)
library(ggiraphExtra)
library(readxl)
# library(sjmisc)
library(plyr)
library(magrittr)
library(stringr)
library(tidyr)
library(radiant)
library(dplyr)
library(lubridate)
library(scales) # to access breaks/formatting functions
library(broom)
# install.packages("stargazer")
library(stargazer)
library(jtools)
# library(lmtest)
library(sandwich)
library(estimatr)
# library(magrittr)
# library(huxtable)
library(latex2exp)
write("", file=".blank")
loadhistory(".blank")
unlink(".blank")
# list.files(path = ".", pattern = NULL, all.files = FALSE,
#            full.names = FALSE, recursive = FALSE,
#            ignore.case = FALSE, include.dirs = FALSE, no.. = FALSE)
# 
# dir(path = ".", pattern = NULL, all.files = FALSE,
#     full.names = FALSE, recursive = FALSE,
#     ignore.case = FALSE, include.dirs = FALSE, no.. = FALSE)
graphics.off()
rm(list = ls())
remove_backticks = function(text){
    text = gsub("([^A-z]+)`", "\\1", text, perl = TRUE)
    text = gsub("`([^A-z]+)", "\\1", text, perl = TRUE)
    text = gsub("(^`)|(`$)", "", text, perl = TRUE)
    return(text)
}
begintable.latex = "
\\begin{table}[H] 
\\centering 
% \\caption{} 
% \\label{} 
\\begin{adjustbox}{center}
\\resizebox{1.0\\linewidth}{!}{"
endtable.latex ="}
\\end{adjustbox}
\\end{table}"
mysignif_tem<-function(x) {
    if ( identical(x, character(0) ) ) {
        return(x)
    } else  {
        x= gsub("," , "", x )
        v = as.numeric(x)
        if ( abs( as.integer(v)) >=100  ) {
            tem = round(v, 0)
        } else {
            tem = signif(v, digits = 3) 
        }
        if ( tem%%1==0 & v-tem!=0 & abs(as.integer(tem)) >= 100 ) { # Case xxx
            # tem = sprintf("%.1f", round(tem, 1))
            tem = sprintf("%1.0f.", round(v, 0))
            tem = sprintf("%1.1f", round(v, 0))
            # } else if ( tem%%10==0 & v-tem!=0 & abs(as.integer(tem)) >= 100 ) { # Case xx0
            #     tem = sprintf("%1.0f.", round(tem, 0))
        } else if ( (tem*100)%%10==0 & v-tem!=0 & abs(as.integer(tem)) < 10 ) { # Case x.x or case x 
            tem = sprintf("%.2f", round(v, 2))
        } else if ( (tem*100)%%1==0 & v-tem!=0 & abs(as.integer(tem)) == 0 ) { # Case 0.xx
            tem = sprintf("%.3f", round(v, 3))
        } else if ( (tem*10)%%10==0 & v-tem!=0 & abs(as.integer(tem)) < 100 ) { # Case xx
            tem = sprintf("%.1f", round(v, 1))
        } 
        return( as.character(tem) ) 
    }
}
mysignif<-function(x) {
    if ( identical(x, character(0) ) ) {
        return(x)
    } else  {
        tem = sapply(x, mysignif_tem, simplify = TRUE, USE.NAMES = FALSE)
    }
    return( tem ) 
}






keepVariables = c("remove_backticks", "data2work", "data2work1", "data2work2", "begintable.latex", "endtable.latex", "keepVariables", "data2work0", "mysignif", "mysignif_tem" ) 
cat("\014")  

# Data loading------------------------------------------------------
data2work0 = readRDS(file = 'data2work20240412_trillion.Rdata')
data2work0 %<>% mutate(top100balance1 = Repo_balances1 + AO_balances1, 
                       logSOFR_IOER_normalized = log(SOFR_IOER - min(data2work0$SOFR_IOER)+ 1) )
data2work1 = 
    apply_labels(
        data2work0,
        IssuanceBil = 'net Bill issuance',
        SecInTotal_bil = 'securities financed in',
        NetPosition_try = 'net Treasuries inventory',
        NetPosition_all = 'net total inventory',
        QuarterEnd = 'quarter-end fixed effect',
        HHI = 'Herfindahl index',
        AfterQE = 'quarter end +1',
        BeforeQE = 'quarter end -1',
        SOFR_IOER = 'SOFR - IOR',
        IssuanceBil_resid = 'Treasury issuance (residual)',
        NetPosition_resid = 'net Treasuries inventory (residual)',
        SecInTotal_bil_resid = 'securities financed in (residual)',
        HHI_resid = 'Herfindahl index (residual)',
        ShortTreasuries_resid = 'Treasury outstanding (residual)',
        ShortTreasuries = 'Treasury outstanding',
        TGCFrepo = 'GCF repo',
        GCF_IOER = 'GCF - IOR',
        total_issuance_bills = 'Bill issuance' ,
        total_issuance_bonds = 'Bond issuance' ,
        total_issuance_notes = 'Note issuance' ,
        total_redemptions_bills = 'Bill redemption' ,
        total_redemptions_bonds = 'Bond redemption' ,
        total_redemptions_notes = 'Note redemption' ,
        total_issuance_all = 'Treasuries issuance' ,
        total_redemptions_all = 'Treasuries redemption' ,
        coupon_issuances = 'Coupon issuance' ,
        logVWATR_IOER = 'log (VWATR - IOR)',
        TGCR_IOER = 'TGCR - IOR',
        logVWATR_IOER_20min = 'log (VWATR_{20min} - IOR)',
        logSOFR_IOER = 'log (SOFR - IOR)',
        Tbills_outstanding = "Tbills outstanding",
        TBills_position = "TBills position",
        Bond_position = "Bonds position",
        Note_position = "Notes position",                            
        TIP_position = "TIPs position",
        AGYMBS_position = "AGYMBS position",
        # 
        deposits = "dealer bank deposits old data",
        DlrDeposits = 	"dealer bank deposits all",
        DlrDepUnins = 	"dealer bank uninsured deposits ",
        OtherDeposits = "other bank deposits all",
        OtherDepUnins = "other bank uninsured deposits ",
        # 
        Repo_balances1 = 'dealer opening balances',
        AO_balances1 = 'other large bank balances',
        Repo_balances2 = 'dealer opening balances',
        AO_balances2 = 'other large bank balances',
        agg_balance_drop1 = 'dealer balance drop',
        agg_balance_drop2 = 'dealer balance drop',
        totalopeningbalances_minus_aggbalancedrop1 = "dealer opening balance - balance drop",
        totalopeningbalances_minus_aggbalancedrop2 = "dealer opening balance - balance drop",
        # 
        tdiff50_recvF1 = "median time of receives from non-dealer banks",
        tdiff25_recvF1 = "Q1 time of receives from non-dealer banks",
        tdiff75_recvF1 = "Q3 time of receives from non-dealer banks",
        tdiff50_all1 =   "median time of receives and sends",
        tdiff25_all1 =   "Q1 time of receives and sends",
        tdiff75_all1 =   "Q3 time of receives and sends",
        tdiff50_send1 =  "median time of sends",
        tdiff25_send1 =  "Q1 time of sends",
        tdiff75_send1 =  "Q3 time of sends",
        tdiff50_recv1 =  "median time of receives",
        tdiff25_recv1 =  "Q1 time of receives",
        tdiff75_recv1 =  "Q3 time of receives",
        tdiff50_recvF2 =  "median time of receives from non-dealer banks",
        tdiff25_recvF2 =  "Q1 time of receives from non-dealer banks",
        tdiff75_recvF2 =  "Q3 time of receives from non-dealer banks",
        tdiff50_all2 =    "median time of receives and sends",
        tdiff25_all2 =    "Q1 time of receives and sends",
        tdiff75_all2 =    "Q3 time of receives and sends",
        tdiff50_send2 =   "median time of sends",
        tdiff25_send2 =   "Q1 time of sends",
        tdiff75_send2 =   "Q3 time of sends",
        tdiff50_recv2 =   "median time of receives",
        tdiff25_recv2 =   "Q1 time of receives",
        tdiff75_recv2 =   "Q3 time of receives", 
        # 
        tdiff25_newIV =   "Q1 time of receives for small banks",
        tdiff50_newIV =   "median time of receives for small banks",
        tdiff75_newIV =   "Q3 time of receives for small banks",
        # 
        Corporation_Income_Taxes_trillion = "corporate tax to US treasury",
        # 
        top100balance1 = "100 large banks",                                   
        logSOFR_IOER_normalized = "log (normalized SOFR-IOR)"
    )

data2work1 %<>% mutate(p25Diff = tdiff25_newIV, p50Diff = tdiff50_newIV, p75Diff = tdiff75_newIV)
cat("\014")

# Summary statistics---------------------------------
# Table 4 and Table 5
rm( list = setdiff( ls(), c(keepVariables, "") ))
cat("\014")  

data2worktem = data2work1 %>%      # we want to show summary statistics in billions. 
    mutate(
        tdiff50_recvF1 =   as.numeric(  tdiff50_recvF1    ),
        tdiff25_recvF1 =   as.numeric(  tdiff25_recvF1    ),
        tdiff75_recvF1 =   as.numeric(  tdiff75_recvF1    ),
        tdiff50_all1   =   as.numeric(  tdiff50_all1      ),
        tdiff25_all1   =   as.numeric(  tdiff25_all1      ),
        tdiff75_all1   =   as.numeric(  tdiff75_all1      ),
        tdiff50_send1  =   as.numeric(  tdiff50_send1     ),
        tdiff25_send1  =   as.numeric(  tdiff25_send1     ),
        tdiff75_send1  =   as.numeric(  tdiff75_send1     ),
        tdiff50_recv1  =   as.numeric(  tdiff50_recv1     ),
        tdiff25_recv1  =   as.numeric(  tdiff25_recv1     ),
        tdiff75_recv1  =   as.numeric(  tdiff75_recv1     ),
        # tdiff50_recvF2 =   as.numeric(  tdiff50_recvF2    ),
        # tdiff25_recvF2 =   as.numeric(  tdiff25_recvF2    ),
        # tdiff75_recvF2 =   as.numeric(  tdiff75_recvF2    ),
        # tdiff50_all2   =   as.numeric(  tdiff50_all2      ),
        # tdiff25_all2   =   as.numeric(  tdiff25_all2      ),
        # tdiff75_all2   =   as.numeric(  tdiff75_all2      ),
        # tdiff50_send2  =   as.numeric(  tdiff50_send2     ),
        # tdiff25_send2  =   as.numeric(  tdiff25_send2     ),
        # tdiff75_send2  =   as.numeric(  tdiff75_send2     ),
        # tdiff50_recv2  =   as.numeric(  tdiff50_recv2     ),
        # tdiff25_recv2  =   as.numeric(  tdiff25_recv2     ),
        # tdiff75_recv2  =   as.numeric(  tdiff75_recv2     ), 
        # 
        Repo_balances1           =   Repo_balances1*1000,
        AO_balances1             =   AO_balances1*1000,
        agg_balance_drop1        =   agg_balance_drop1*1000,
        #                       
        Tbills_outstanding       =   Tbills_outstanding*1000,
        #                          
        total_issuance_bills     =   total_issuance_bills*1000,
        coupon_issuances         =   coupon_issuances*1000,
        TBills_position          =   TBills_position*1000,
        total_redemptions_all    =   total_redemptions_all*1000,
        NetPosition_try          =   NetPosition_try*1000,
        #                          
        #                          
        total_issuance_notes     =   total_issuance_notes*1000,
        total_issuance_bonds     =   total_issuance_bonds*1000,
        total_issuance_all       =   total_issuance_all*1000,
        #                                         
        #                                          
        total_redemptions_bills  =   total_redemptions_bills*1000,
        total_redemptions_bonds  =   total_redemptions_bonds*1000,
        total_redemptions_notes  =   total_redemptions_notes*1000,
        coupon_redemptions       =   coupon_redemptions*1000,
        #                                             
        #                                            
        Bond_position            =   Bond_position*1000,
        Note_position            =   Note_position*1000,
        TIP_position             =   TIP_position*1000,
        AGYMBS_position          =   AGYMBS_position*1000,
        NetPosition_all          =   NetPosition_all*1000,
        Corporation_Income_Taxes_trillion = Corporation_Income_Taxes_trillion*1000,
        top100balance1 = top100balance1*1000,
        DlrDepUnins = DlrDepUnins * 1000
    )


data2worktem1 = data2worktem %>% 
    select(Date,
           Repo_balances1  ,
           AO_balances1  ,
           top100balance1,
           # agg_balance_drop1  ,
           #
           Tbills_outstanding  ,
           #
           total_issuance_bills  ,
           coupon_issuances  ,
           TBills_position  ,
           total_redemptions_all  ,
           NetPosition_try  ,
           tdiff50_recv1,
           tdiff50_send1,
           tdiff25_newIV,
           tdiff50_newIV  ,
           tdiff75_newIV    ,
           SOFR_IOER, 
           # GCF_IOER,
           # TGCR_IOER,  
           total_issuance_all, 
           QuarterEnd ,
           Corporation_Income_Taxes_trillion,
           DlrDepUnins,
           logSOFR_IOER_normalized
    ) %>% mutate(#Ave_tdiff50_recv1_3 = zoo::rollmean(tdiff50_recv1, k = 3, align = "right", fill = NA, na.rm = TRUE),
                 # Ave_tdiff50_recv1_5 = zoo::rollmean(tdiff50_recv1, k = 5, align = "right", fill = NA, na.rm = TRUE),
                 # Ave_tdiff50_recv1_7 = zoo::rollmean(tdiff50_recv1, k = 7, align = "right", fill = NA, na.rm = TRUE),
                 Ave_tdiff50_recv1_9 = zoo::rollmean(tdiff50_recv1, k = 9, align = "right", fill = NA, na.rm = TRUE),
                 # Ave_tdiff50_recv1_11 = zoo::rollmean(tdiff50_recv1, k = 11, align = "right", fill = NA, na.rm = TRUE),
                 # Ave_tdiff50_recv1_13 = zoo::rollmean(tdiff50_recv1, k = 13, align = "right", fill = NA, na.rm = TRUE),
                 # Ave_tdiff50_recv1_15 = zoo::rollmean(tdiff50_recv1, k = 15, align = "right", fill = NA, na.rm = TRUE),
                 # Ave_tdiff50_recv1_17 = zoo::rollmean(tdiff50_recv1, k = 17, align = "right", fill = NA, na.rm = TRUE),
                 # Ave_tdiff50_recv1_19 = zoo::rollmean(tdiff50_recv1, k = 19, align = "right", fill = NA, na.rm = TRUE),
                 # Ave_tdiff50_recv1_21 = zoo::rollmean(tdiff50_recv1, k = 21, align = "right", fill = NA, na.rm = TRUE),
                 # AveBalancesPreviousWeeklag3 = lag(Ave_tdiff50_recv1_3, 2),
                 # AveBalancesPreviousWeeklag5 = lag(Ave_tdiff50_recv1_5, 2),
                 # AveBalancesPreviousWeeklag7 = lag(Ave_tdiff50_recv1_7, 2),
                 AveBalancesPreviousWeeklag9 = lag(Ave_tdiff50_recv1_9, 2)
                 # AveBalancesPreviousWeeklag11 = lag(Ave_tdiff50_recv1_11, 2),
                 # AveBalancesPreviousWeeklag13 = lag(Ave_tdiff50_recv1_13, 2),
                 # AveBalancesPreviousWeeklag15 = lag(Ave_tdiff50_recv1_15, 2),
                 # AveBalancesPreviousWeeklag17 = lag(Ave_tdiff50_recv1_17, 2),
                 # AveBalancesPreviousWeeklag19 = lag(Ave_tdiff50_recv1_19, 2),
                 # AveBalancesPreviousWeeklag21 = lag(Ave_tdiff50_recv1_21, 2),
                 # lag2_Repo_balances1 = lag(Repo_balances1, 2)
    )

data2worktem1 = 
    apply_labels(
        data2worktem1,
        Repo_balances1	=	paste0(		var_lab(data2worktem1$Repo_balances1)		, "	($ billions)"	),
        AO_balances1	=	paste0(		var_lab(data2worktem1$AO_balances1)		, "	($ billions)"	),
        top100balance1	=	paste0(		var_lab(data2worktem1$top100balance1)		, "	($ billions)"	),
        Tbills_outstanding	=	paste0(		var_lab(data2worktem1$Tbills_outstanding)		, "	($ billions)"	),
        #	=	paste0(		var_lab(data2worktem1$#)		, "	($ billions)"	),
        total_issuance_bills	=	paste0(		var_lab(data2worktem1$total_issuance_bills)		, "	($ billions)"	),
        coupon_issuances	=	paste0(		var_lab(data2worktem1$coupon_issuances)		, "	($ billions)"	),
        TBills_position	=	paste0(		var_lab(data2worktem1$TBills_position)		, "	($ billions)"	),
        total_redemptions_all	=	paste0(		var_lab(data2worktem1$total_redemptions_all)		, "	($ billions)"	),
        NetPosition_try	=	paste0(		var_lab(data2worktem1$NetPosition_try)		, "	($ billions)"	),
        tdiff50_recv1	=	paste0(		var_lab(data2worktem1$tdiff50_recv1)		, "	(minutes)"	),
        tdiff50_send1	=	paste0(		var_lab(data2worktem1$tdiff50_send1)		, "	(minutes)"	),
        # 
        tdiff25_newIV	=	paste0(		var_lab(data2worktem1$tdiff25_newIV)		, "	(minutes)"	),
        tdiff50_newIV	=	paste0(		var_lab(data2worktem1$tdiff50_newIV)		, "	(minutes)"	),
        tdiff75_newIV	=	paste0(		var_lab(data2worktem1$tdiff75_newIV)		, "	(minutes)"	),
        # 
        SOFR_IOER	=	paste0(		var_lab(data2worktem1$SOFR_IOER)		, "	(basis points)"	),
        # GCF_IOER	=	paste0(		var_lab(data2worktem1$GCF_IOER)		, "	(basis points)"	),
        #TGCR_IOER	=	paste0(		var_lab(data2worktem1$#TGCR_IOER)		, "	(basis points)"	),
        total_issuance_all	=	paste0(		var_lab(data2worktem1$total_issuance_all)		, "	($ billions)"	),
        QuarterEnd	=	paste0(		var_lab(data2worktem1$QuarterEnd)		, "	"	),
        Corporation_Income_Taxes_trillion	=	paste0(		var_lab(data2worktem1$Corporation_Income_Taxes_trillion)		, "	($ billions)"	),
        DlrDepUnins	=	paste0(		var_lab(data2worktem1$DlrDepUnins)		, "	($ billions)"	),
        logSOFR_IOER_normalized	=	paste0(		var_lab(data2worktem1$logSOFR_IOER_normalized)		, "	"	),
        AveBalancesPreviousWeeklag9 = "Average Payment Delay t-10 to t-2"
    )

data2work7 = names2labels(data2worktem1)

data2work7 %<>% select(-Date)

res = remove_backticks( 
    stargazer(data2work7,
              digits = 1, summary.stat = c("n",	"mean",	"sd",	"min",	"p25",	"median",	"p75",	"max"                )

    ))
res[grepl("\\label\\{\\}",res)] = "%"
res[grepl("\\caption\\{",res)] = "%"

res[grepl("Table created by stargazer v.5.2.2 by Marek Hlavac",res)] = ""
res[grepl("\\% Date and time:",res)] = ""
res[grepl("\\% Requires LaTeX packages: dcolumn",res)] = ""
# 
res[grepl("begin\\{table\\}",res)] = begintable.latex
res[grepl("end\\{table\\}",res)] = endtable.latex
cat("\014")  

write_clip(res)
writeLines(res)
cat("\014")  
# end


# **********************************************************
rm( list = setdiff( ls(), c(keepVariables, "data2worktem") ))

data2worktem1 = data2worktem %>% 
    select(#SOFR	,
           # TGCR	,
           # GCF	,
           IOER	,
           # BGCR	,
           TGCR_IOER	,
           GCF_IOER	,
           # BGCR_IOER	,
           total_issuance_bonds	,
           total_issuance_notes	,
           total_redemptions_bills	,
           total_redemptions_bonds	,
           total_redemptions_notes	,
           # coupon_redemptions	,
           # Bond_position	,
           # Note_position	,
           # TIP_position	,
           # AGYMBS_position	,
           NetPosition_all	,
           # BeforeQE	,
           # AfterQE	,
           # agg_balance_drop1	,
           # tdiff50_recvF1	,
           # tdiff25_recvF1	,
           # tdiff75_recvF1	,
           tdiff50_all1	,
           tdiff25_all1	,
           tdiff75_all1	,
           tdiff25_send1	,
           tdiff75_send1	,
           tdiff25_recv1	,
           tdiff75_recv1	
    )
# write_clip(colnames(data2worktem1))

data2worktem1 = 
    apply_labels(
        data2worktem1,
        # SOFR	=	paste0(		var_lab(data2worktem1$SOFR)		, "SOFR	(basis points)"	),
        # TGCR	=	paste0(		var_lab(data2worktem1$# TGCR)		, "	(basis points)"	),
        # GCF	=	paste0(		var_lab(data2worktem1$# GCF)		, "	(basis points)"	),
        IOER	=	paste0(		var_lab(data2worktem1$IOER)		, "IOR	(basis points)"	),
        # BGCR	=	paste0(		var_lab(data2worktem1$# BGCR)		, "	(basis points)"	),
        TGCR_IOER	=	paste0(		var_lab(data2worktem1$TGCR_IOER)		, "	(basis points)"	),
        GCF_IOER	=	paste0(		var_lab(data2worktem1$GCF_IOER)		, "	(basis points)"	),
        # BGCR_IOER	=	paste0(		var_lab(data2worktem1$# BGCR_IOER)		, "	(basis points)"	),
        total_issuance_bonds	=	paste0(		var_lab(data2worktem1$total_issuance_bonds)		, "	($ billions)"	),
        total_issuance_notes	=	paste0(		var_lab(data2worktem1$total_issuance_notes)		, "	($ billions)"	),
        total_redemptions_bills	=	paste0(		var_lab(data2worktem1$total_redemptions_bills)		, "	($ billions)"	),
        total_redemptions_bonds	=	paste0(		var_lab(data2worktem1$total_redemptions_bonds)		, "	($ billions)"	),
        total_redemptions_notes	=	paste0(		var_lab(data2worktem1$total_redemptions_notes)		, "	($ billions)"	),
        # coupon_redemptions	=	paste0(		var_lab(data2worktem1$coupon_redemptions)		, "	($ billions)"	),
        # Bond_position	=	paste0(		var_lab(data2worktem1$# Bond_position)		, "	($ billions)"	),
        # Note_position	=	paste0(		var_lab(data2worktem1$# Note_position)		, "	($ billions)"	),
        # TIP_position	=	paste0(		var_lab(data2worktem1$# TIP_position)		, "	($ billions)"	),
        # AGYMBS_position	=	paste0(		var_lab(data2worktem1$# AGYMBS_position)		, "	($ billions)"	),
        NetPosition_all	=	paste0(		var_lab(data2worktem1$NetPosition_all)		, "	($ billions)"	),
        # BeforeQE	=	paste0(		var_lab(data2worktem1$# BeforeQE)		, "	"	),
        # AfterQE	=	paste0(		var_lab(data2worktem1$# AfterQE)		, "	"	),
        # agg_balance_drop1	=	paste0(		var_lab(data2worktem1$# agg_balance_drop1)		, "	($ billions)"	),
        # tdiff50_recvF1	=	paste0(		var_lab(data2worktem1$# tdiff50_recvF1)		, "	"	),
        # tdiff25_recvF1	=	paste0(		var_lab(data2worktem1$# tdiff25_recvF1)		, "	(minutes)"	),
        # tdiff75_recvF1	=	paste0(		var_lab(data2worktem1$# tdiff75_recvF1)		, "	(minutes)"	),
        tdiff50_all1	=	paste0(		var_lab(data2worktem1$tdiff50_all1)		, "	(minutes)"	),
        tdiff25_all1	=	paste0(		var_lab(data2worktem1$tdiff25_all1)		, "	(minutes)"	),
        tdiff75_all1	=	paste0(		var_lab(data2worktem1$tdiff75_all1)		, "	(minutes)"	),
        tdiff25_send1	=	paste0(		var_lab(data2worktem1$tdiff25_send1)		, "	(minutes)"	),
        tdiff75_send1	=	paste0(		var_lab(data2worktem1$tdiff75_send1)		, "	(minutes)"	),
        tdiff25_recv1	=	paste0(		var_lab(data2worktem1$tdiff25_recv1)		, "	(minutes)"	),
        tdiff75_recv1	=	paste0(		var_lab(data2worktem1$tdiff75_recv1)		, "	(minutes)"	)
    )

data2work7 = names2labels(data2worktem1)


res = remove_backticks( 
    stargazer(data2work7,
              digits = 1
    ))
res[grepl("\\label\\{\\}",res)] = "%"
res[grepl("\\caption\\{",res)] = "%"

res[grepl("Table created by stargazer v.5.2.2 by Marek Hlavac",res)] = ""
res[grepl("\\% Date and time:",res)] = ""
res[grepl("\\% Requires LaTeX packages: dcolumn",res)] = ""
# regmatches(res, gregexpr("(-?)(\\d+(\\.\\d+)?)", res)) =  lapply(regmatches(res, gregexpr("(-?)(\\d+(\\.\\d+)?)", res) ) , mysignif )
# 
# res[grepl("begin\\{table\\}",res)] = "~\\newpage
# \\begin{table}[H] 
# \\centering 
# % \\caption{} 
# % \\label{} 
# \\begin{adjustbox}{center}
# \\resizebox{0.8\\linewidth}{!}{"
# 
res[grepl("begin\\{table\\}",res)] = begintable.latex
res[grepl("end\\{table\\}",res)] = endtable.latex
cat("\014")  

write_clip(res)
writeLines(res)
cat("\014")  














































# ----------------------Regression-----------------------------
rm( list = setdiff( ls(), c(keepVariables, "") ))

cat("\014")  
# ---------------------Table 6-----------------------------------------------------------------------
# First definition
rm( list = setdiff( ls(), c(keepVariables, "") ))
# data2work = data2work1 %>% mutate(top100balance1 = Repo_balances1 + AO_balances1, 
#                                   logSOFR_IOER_normalized = log(SOFR_IOER - min(data2work1$SOFR_IOER)+ 1) )
# data2work = apply_labels(data2work, top100balance1 = "100 large banks",                                   
#                           logSOFR_IOER_normalized = "log (normalized SOFR-IOR)")
data2work = data2work1 
regressor_list = c(
    'tdiff50_recv1 ~ top100balance1', 
    'tdiff50_recv1 ~ Repo_balances1 + AO_balances1',
    'tdiff50_recv1 ~ Repo_balances1 + AO_balances1 + logSOFR_IOER_normalized + total_issuance_all + total_redemptions_all + DlrDepUnins',
    'tdiff50_recv1 ~ Repo_balances1 + AO_balances1 + SOFR_IOER + total_issuance_all + total_redemptions_all + DlrDepUnins',
    
    'tdiff50_send1 ~ top100balance1', 
    'tdiff50_send1 ~ Repo_balances1 + AO_balances1',
    'tdiff50_send1 ~ Repo_balances1 + AO_balances1 + logSOFR_IOER_normalized + total_issuance_all + total_redemptions_all + DlrDepUnins',
    'tdiff50_send1 ~ Repo_balances1 + AO_balances1 + SOFR_IOER + total_issuance_all + total_redemptions_all + DlrDepUnins',
    'tdiff50_send1 ~ Repo_balances1 + AO_balances1 + SOFR_IOER + total_issuance_all + total_redemptions_all + tdiff50_recv1 + DlrDepUnins'
    )

dep_list = c('tdiff50_recv1', 'tdiff50_send1')

reg_list = vector(mode = 'list', length = length(regressor_list) * 1 )
robust_se = vector(mode = 'list', length = length(regressor_list) * 1  )
x_vars = c()
i = 1

for (a1 in regressor_list) {
    tmp_formula_str = a1
    commandline = paste0('use_labels(data2work, lm(formula(', tmp_formula_str, '), data = data2work ) )' )
    tmp_model = eval(parse(text = commandline))
    # tmp_model = lm(as.formula(tmp_formula_str), data = data2work )
    # tmp_model = lm(as.formula(tmp_formula_str), data = Data2regress)
    cov1         = vcovHC(tmp_model, type = "HC1")
    robust_se[[i]] = sqrt(diag(cov1))
    reg_list[[i]] = tmp_model
    a2 = if_else(i<=4, dep_list[1], dep_list[2])
    x_vars = c(x_vars, a2)
    i = i + 1
}





note.latex = paste0(
    "\\multicolumn{", 
    length(reg_list)+1,
    "}{l}{{Note:} Standard errors are adjusted for heteroskedasticity. $^{*}p<0.1$; $^{**}p<0.05$; $^{***}p<0.01$.} \\\\ \\multicolumn{",
    length(reg_list)+1,
    "}{l}{ \\qquad \\quad Constant included for each specification.}")

res = remove_backticks( 
    stargazer(
        reg_list, 
        title="Results", 
        align=TRUE, 
        omit.stat = c("f"), 
        dep.var.caption  =  paste0("Dependent variable: ", eval(parse(text = paste0('var_lab(data2work$',dep_list[1],')'))) ),
        dep.var.labels.include = FALSE,
        df = FALSE,
        digits = 5,
        # notes        = TeX("Sometimes you just have to start over."), 
        # notes.append = FALSE,
        se = robust_se))
res[grepl("Note",res)] = note.latex
res[grepl("\\label\\{\\}",res)] = "%"
res[grepl("\\caption\\{",res)] = "%"

res[grepl("Table created by stargazer",res)] = ""
res[grepl("\\% Date and time:",res)] = ""
res[grepl("\\% Requires LaTeX packages: dcolumn",res)] = ""

regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res)) =  lapply(regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res) ) , mysignif )
# 
# res[grepl("begin\\{table\\}",res)] = "~\\newpage
# \\begin{table}[H] 
# \\centering 
# % \\caption{} 
# % \\label{} 
# \\begin{adjustbox}{center}
# \\resizebox{0.8\\linewidth}{!}{"
# 
res[grepl("begin\\{table\\}",res)] = begintable.latex
res[grepl("end\\{table\\}",res)] = endtable.latex
res = gsub("R\\$\\^\\{2\\}\\$", "\\$R\\^\\{2\\}\\$", res)
cat("\014")  
write_clip(res)
writeLines(res)


# ---------------------Main Regression -----------------------------------------------------------------------
# First definition
# Table 1

rm( list = setdiff( ls(), c(keepVariables, "") ))
data2work = data2work1

regressor_list = c(
    'Repo_balances1', # 'tdiff50_recv1',
    'Repo_balances1 + QuarterEnd',
    'tdiff50_recv1 + QuarterEnd',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding',
    # 'Repo_balances1 + QuarterEnd + Tbills_outstanding',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + DlrDepUnins',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + DlrDepUnins',
    # 'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins + Corporation_Income_Taxes_trillion'
)

dep_list = c('SOFR_IOER')

reg_list = vector(mode = 'list', length = length(regressor_list) * length(dep_list) )
robust_se = vector(mode = 'list', length = length(regressor_list) * length(dep_list)  )
x_vars = c()
i = 1

for (a1 in regressor_list) {
    for (a2 in dep_list) {
        tmp_formula_str = paste0(a2, ' ~ ', a1)
        commandline = paste0('use_labels(data2work, lm(formula(', tmp_formula_str, '), data = data2work ) )' )
        tmp_model = eval(parse(text = commandline))
        # tmp_model = lm(as.formula(tmp_formula_str), data = data2work )
        # tmp_model = lm(as.formula(tmp_formula_str), data = Data2regress)
        cov1         = vcovHC(tmp_model, type = "HC1")
        robust_se[[i]] = sqrt(diag(cov1))
        reg_list[[i]] = tmp_model
        x_vars = c(x_vars, a2)
        i = i + 1
    }
}
note.latex = paste0(
    "\\multicolumn{", 
    length(reg_list)+1,
    "}{l}{{Note:} Standard errors are adjusted for heteroskedasticity. $^{*}p<0.1$; $^{**}p<0.05$; $^{***}p<0.01$.} \\\\ \\multicolumn{",
    length(reg_list)+1,
    "}{l}{ \\qquad \\quad Constant included for each specification.}")


res = remove_backticks( 
    stargazer(
        reg_list, 
        title="Results", 
        align=TRUE, 
        omit.stat = c("f"), 
        dep.var.caption  =  paste0("Dependent variable: ", eval(parse(text = paste0('var_lab(data2work$',dep_list[1],')'))) ),
        dep.var.labels.include = FALSE,
        df = FALSE,
        digits = 5,
        # notes        = TeX("Sometimes you just have to start over."), 
        # notes.append = FALSE,
        se = robust_se))
res[grepl("Note",res)] = note.latex
res[grepl("\\label\\{\\}",res)] = "%"
res[grepl("\\caption\\{",res)] = "%"

res[grepl("Table created by stargazer",res)] = ""
res[grepl("\\% Date and time:",res)] = ""
res[grepl("\\% Requires LaTeX packages: dcolumn",res)] = ""

regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res)) =  lapply(regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res) ) , mysignif )
# 
res[grepl("begin\\{table\\}",res)] = begintable.latex
res[grepl("end\\{table\\}",res)] = endtable.latex
res = gsub("R\\$\\^\\{2\\}\\$", "\\$R\\^\\{2\\}\\$", res)
cat("\014")  
write_clip(res)
writeLines(res)
















# ---------------------Main Regression Shorter Sample 1-----------------------------------------------------------------------
# First definition
rm( list = setdiff( ls(), c(keepVariables, "") ))
data2work = data2work1 %>% filter(Date <= ymd('2020-10-30'))
# colnames(data2work)

regressor_list = c(
    'Repo_balances1', # 'tdiff50_recv1',
    'Repo_balances1 + QuarterEnd',
    'tdiff50_recv1 + QuarterEnd',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding',
    # 'Repo_balances1 + QuarterEnd + Tbills_outstanding',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + DlrDepUnins',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + DlrDepUnins',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins + Corporation_Income_Taxes_trillion'
)

dep_list = c('SOFR_IOER')

reg_list = vector(mode = 'list', length = length(regressor_list) * length(dep_list) )
robust_se = vector(mode = 'list', length = length(regressor_list) * length(dep_list)  )
x_vars = c()
i = 1

for (a1 in regressor_list) {
    for (a2 in dep_list) {
        tmp_formula_str = paste0(a2, ' ~ ', a1)
        commandline = paste0('use_labels(data2work, lm(formula(', tmp_formula_str, '), data = data2work ) )' )
        tmp_model = eval(parse(text = commandline))
        # tmp_model = lm(as.formula(tmp_formula_str), data = data2work )
        # tmp_model = lm(as.formula(tmp_formula_str), data = Data2regress)
        cov1         = vcovHC(tmp_model, type = "HC1")
        robust_se[[i]] = sqrt(diag(cov1))
        reg_list[[i]] = tmp_model
        x_vars = c(x_vars, a2)
        i = i + 1
    }
}
note.latex = paste0(
    "\\multicolumn{", 
    length(reg_list)+1,
    "}{l}{{Note:} Standard errors are adjusted for heteroskedasticity. $^{*}p<0.1$; $^{**}p<0.05$; $^{***}p<0.01$.} \\\\ \\multicolumn{",
    length(reg_list)+1,
    "}{l}{ \\qquad \\quad Constant included for each specification.}")

res = remove_backticks( 
    stargazer(
        reg_list, 
        title="Results", 
        align=TRUE, 
        omit.stat = c("f"), 
        dep.var.caption  =  paste0("Dependent variable: ", eval(parse(text = paste0('var_lab(data2work$',dep_list[1],')'))) ),
        dep.var.labels.include = FALSE,
        df = FALSE,
        digits = 5,
        # notes        = TeX("Sometimes you just have to start over."), 
        # notes.append = FALSE,
        se = robust_se))
res[grepl("Note",res)] = note.latex
res[grepl("\\label\\{\\}",res)] = "%"
res[grepl("\\caption\\{",res)] = "%"

res[grepl("Table created by stargazer",res)] = ""
res[grepl("\\% Date and time:",res)] = ""
res[grepl("\\% Requires LaTeX packages: dcolumn",res)] = ""

regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res)) =  lapply(regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res) ) , mysignif )
# 
res[grepl("begin\\{table\\}",res)] = begintable.latex
res[grepl("end\\{table\\}",res)] = endtable.latex
res = gsub("R\\$\\^\\{2\\}\\$", "\\$R\\^\\{2\\}\\$", res)
cat("\014")  
write_clip(res)
writeLines(res)




# ---------------------Main Regression Shorter Sample 2-----------------------------------------------------------------------
# First definition
rm( list = setdiff( ls(), c(keepVariables, "") ))
data2work = data2work1 %>% filter(Date <= ymd('2020-3-1'))
# colnames(data2work)

regressor_list = c(
    'Repo_balances1', # 'tdiff50_recv1',
    'Repo_balances1 + QuarterEnd',
    'tdiff50_recv1 + QuarterEnd',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding',
    # 'Repo_balances1 + QuarterEnd + Tbills_outstanding',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + DlrDepUnins',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + DlrDepUnins',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins + Corporation_Income_Taxes_trillion'
)
# 

dep_list = c('SOFR_IOER')

reg_list = vector(mode = 'list', length = length(regressor_list) * length(dep_list) )
robust_se = vector(mode = 'list', length = length(regressor_list) * length(dep_list)  )
x_vars = c()
i = 1

for (a1 in regressor_list) {
    for (a2 in dep_list) {
        tmp_formula_str = paste0(a2, ' ~ ', a1)
        commandline = paste0('use_labels(data2work, lm(formula(', tmp_formula_str, '), data = data2work ) )' )
        tmp_model = eval(parse(text = commandline))
        # tmp_model = lm(as.formula(tmp_formula_str), data = data2work )
        # tmp_model = lm(as.formula(tmp_formula_str), data = Data2regress)
        cov1         = vcovHC(tmp_model, type = "HC1")
        robust_se[[i]] = sqrt(diag(cov1))
        reg_list[[i]] = tmp_model
        x_vars = c(x_vars, a2)
        i = i + 1
    }
}
note.latex = paste0(
    "\\multicolumn{", 
    length(reg_list)+1,
    "}{l}{{Note:} Standard errors are adjusted for heteroskedasticity. $^{*}p<0.1$; $^{**}p<0.05$; $^{***}p<0.01$.} \\\\ \\multicolumn{",
    length(reg_list)+1,
    "}{l}{ \\qquad \\quad Constant included for each specification.}")

res = remove_backticks( 
    stargazer(
        reg_list, 
        title="Results", 
        align=TRUE, 
        omit.stat = c("f"), 
        dep.var.caption  =  paste0("Dependent variable: ", eval(parse(text = paste0('var_lab(data2work$',dep_list[1],')'))) ),
        dep.var.labels.include = FALSE,
        df = FALSE,
        digits = 5,
        # notes        = TeX("Sometimes you just have to start over."), 
        # notes.append = FALSE,
        se = robust_se))
res[grepl("Note",res)] = note.latex
res[grepl("\\label\\{\\}",res)] = "%"
res[grepl("\\caption\\{",res)] = "%"

res[grepl("Table created by stargazer",res)] = ""
res[grepl("\\% Date and time:",res)] = ""
res[grepl("\\% Requires LaTeX packages: dcolumn",res)] = ""

regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res)) =  lapply(regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res) ) , mysignif )
# 
# 
res[grepl("begin\\{table\\}",res)] = begintable.latex
res[grepl("end\\{table\\}",res)] = endtable.latex
res = gsub("R\\$\\^\\{2\\}\\$", "\\$R\\^\\{2\\}\\$", res)
cat("\014")  
write_clip(res)
writeLines(res)




# ---------------------Main Regression without spike days-----------------------------------------------------------------------
# First definition
rm( list = setdiff( ls(), c(keepVariables, "") ))
cat("\014")  

tem1 = zoo::na.locf(data2work1, na.rm = FALSE)
tem2 = tem1 %>% 
    mutate(SOFR_IOER = pmin(SOFR_IOER, 30),
           GCF_IOER  = pmin(GCF_IOER,  30),
           TGCR_IOER = pmin(TGCR_IOER, 30)   ) %>%
    mutate(smooth_SOFR_IOER =  zoo::rollmean( lag(SOFR_IOER), 15, fill=NA, align="right"),
           smooth_GCF_IOER  =  zoo::rollmean( lag(GCF_IOER), 15, fill=NA, align="right"),
           smooth_TGCR_IOER =  zoo::rollmean( lag(TGCR_IOER), 15, fill=NA, align="right"), 
           dif_SOFR_IOER = SOFR_IOER - smooth_SOFR_IOER, 
           dif_GCF_IOER  = GCF_IOER  - smooth_GCF_IOER ,
           dif_TGCR_IOER = TGCR_IOER - smooth_TGCR_IOER) %>%
    select(Date, smooth_SOFR_IOER, smooth_GCF_IOER, smooth_TGCR_IOER, dif_SOFR_IOER, dif_GCF_IOER, dif_TGCR_IOER)

data2work = data2work1 %>% 
    left_join(tem2, by = 'Date') %>%
    mutate(critical_flag1 = ifelse( (dif_SOFR_IOER >= 15 | dif_GCF_IOER >= 15 | dif_TGCR_IOER >= 15 | SOFR_IOER >=20 | GCF_IOER >= 30 | TGCR_IOER >= 20) & !is.na(SOFR_IOER), 1, 0)  ) %>%
    mutate(critical_flag = ifelse(is.na(critical_flag1), 0, critical_flag1) ) %>% select(-critical_flag1)
cat("\014")

data2work %<>% filter(critical_flag != 1 )

regressor_list = c(
    'Repo_balances1', # 'tdiff50_recv1',
    'Repo_balances1 + QuarterEnd',
    'tdiff50_recv1 + QuarterEnd',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding',
    # 'Repo_balances1 + QuarterEnd + Tbills_outstanding',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + DlrDepUnins',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + DlrDepUnins',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins + Corporation_Income_Taxes_trillion'
)

dep_list = c('SOFR_IOER')

reg_list = vector(mode = 'list', length = length(regressor_list) * length(dep_list) )
robust_se = vector(mode = 'list', length = length(regressor_list) * length(dep_list)  )
x_vars = c()
i = 1

for (a1 in regressor_list) {
    for (a2 in dep_list) {
        tmp_formula_str = paste0(a2, ' ~ ', a1)
        commandline = paste0('use_labels(data2work, lm(formula(', tmp_formula_str, '), data = data2work ) )' )
        tmp_model = eval(parse(text = commandline))
        cov1         = vcovHC(tmp_model, type = "HC1")
        robust_se[[i]] = sqrt(diag(cov1))
        reg_list[[i]] = tmp_model
        x_vars = c(x_vars, a2)
        i = i + 1
    }
}
note.latex = paste0(
    "\\multicolumn{", 
    length(reg_list)+1,
    "}{l}{{Note:} Standard errors are adjusted for heteroskedasticity. $^{*}p<0.1$; $^{**}p<0.05$; $^{***}p<0.01$.} \\\\ \\multicolumn{",
    length(reg_list)+1,
    "}{l}{ \\qquad \\quad Constant included for each specification.}")

res = remove_backticks( 
    stargazer(
        reg_list, 
        title="Results", 
        align=TRUE, 
        omit.stat = c("f"), 
        dep.var.caption  =  paste0("Dependent variable: ", eval(parse(text = paste0('var_lab(data2work$',dep_list[1],')'))) ),
        dep.var.labels.include = FALSE,
        df = FALSE,
        digits = 5,
        # notes        = TeX("Sometimes you just have to start over."), 
        # notes.append = FALSE,
        se = robust_se))
res[grepl("Note",res)] = note.latex
res[grepl("\\label\\{\\}",res)] = "%"
res[grepl("\\caption\\{",res)] = "%"

res[grepl("Table created by stargazer",res)] = ""
res[grepl("\\% Date and time:",res)] = ""
res[grepl("\\% Requires LaTeX packages: dcolumn",res)] = ""

regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res)) =  lapply(regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res) ) , mysignif )
# 
# 
res[grepl("begin\\{table\\}",res)] = begintable.latex
res[grepl("end\\{table\\}",res)] = endtable.latex
res = gsub("R\\$\\^\\{2\\}\\$", "\\$R\\^\\{2\\}\\$", res)
cat("\014")  
write_clip(res)
writeLines(res)








# -----------------Probit Regression -----------------------------------------
# Table 10 and Table 11
rm( list = setdiff( ls(), c(keepVariables, "") ))
cat("\014")  

tem1 = zoo::na.locf(data2work1, na.rm = FALSE)
tem2 = tem1 %>% 
    mutate(SOFR_IOER2 = SOFR_IOER,
           GCF_IOER2 = GCF_IOER,
           TGCR_IOER2 = TGCR_IOER,
           SOFR_IOER = pmin(SOFR_IOER, 30),
           GCF_IOER  = pmin(GCF_IOER,  30),
           TGCR_IOER = pmin(TGCR_IOER, 30)
           ) %>%
    mutate(smooth_SOFR_IOER =  zoo::rollmean( lag(SOFR_IOER), 15, fill=NA, align="right"),
           smooth_GCF_IOER  =  zoo::rollmean( lag(GCF_IOER), 15, fill=NA, align="right"),
           smooth_TGCR_IOER =  zoo::rollmean( lag(TGCR_IOER), 15, fill=NA, align="right"), 
           dif_SOFR_IOER = SOFR_IOER2 - smooth_SOFR_IOER, 
           dif_GCF_IOER  = GCF_IOER2  - smooth_GCF_IOER ,
           dif_TGCR_IOER = TGCR_IOER2 - smooth_TGCR_IOER) %>%
    select(Date, smooth_SOFR_IOER, smooth_GCF_IOER, smooth_TGCR_IOER, dif_SOFR_IOER, dif_GCF_IOER, dif_TGCR_IOER)

data2work = data2work1 %>% 
    left_join(tem2, by = 'Date') %>%
    mutate(critical_flag1 = ifelse( (dif_SOFR_IOER >= 15 | dif_GCF_IOER >= 15 | dif_TGCR_IOER >= 15 | SOFR_IOER >=20 | GCF_IOER >= 30 | TGCR_IOER >= 20) & !is.na(SOFR_IOER), 1, 0)  ) %>%
    mutate(critical_flag = ifelse(is.na(critical_flag1), 0, critical_flag1) ) %>% select(-critical_flag1)
# mutate(critical_flag = if_else( (SOFR_IOER >25 | GCF_IOER >25 ), 1, 0 ))
cat("\014")


# colnames(data2worktem1)
data2worktem1 = data2work %>% filter(critical_flag==1) %>%
    select(Date, 
           SOFR_IOER, 
           GCF_IOER,
           TGCR_IOER,
           Repo_balances1,
           AO_balances1,
           total_issuance_all,
           total_redemptions_all
           ) %>% 
    # Change units to billion dollars
    mutate(Repo_balances1	=	Repo_balances1	*	1000	,
           AO_balances1	=	AO_balances1	*	1000	,
           total_issuance_all	=	total_issuance_all	*	1000	,
           total_redemptions_all	=	total_redemptions_all	*	1000	)



# Table 10
data2work7.2 = names2labels(data2worktem1)

res = remove_backticks( 
    stargazer(data2work7.2, summary = FALSE,
              digits = 2))

res[grepl("\\label\\{\\}",res)] = "%"
res[grepl("\\caption\\{",res)] = "%"

res[grepl("Table created by stargazer v.5.2.2 by Marek Hlavac",res)] = ""
res[grepl("\\% Date and time:",res)] = ""
res[grepl("\\% Requires LaTeX packages: dcolumn",res)] = ""
# 
res[grepl("begin\\{table\\}",res)] = begintable.latex
res[grepl("end\\{table\\}",res)] = endtable.latex
cat("\014")  

write_clip(res)
writeLines(res)
cat("\014")  






# Table 11
regressor_list = c(
    'tdiff50_recv1', 'Repo_balances1',
    'tdiff50_recv1 + Repo_balances1',    'coupon_issuances', 
    'tdiff50_recv1 + coupon_issuances',
    'tdiff50_recv1 + Repo_balances1 + coupon_issuances',
    'tdiff50_recv1 + Repo_balances1 + coupon_issuances + QuarterEnd'
)

dep_list = c('critical_flag')

reg_list = vector(mode = 'list', length = length(regressor_list) * length(dep_list) )
robust_se = vector(mode = 'list', length = length(regressor_list) * length(dep_list)  )
x_vars = c()
i = 1

for (a1 in regressor_list) {
    for (a2 in dep_list) {
        tmp_formula_str = paste0(a2, ' ~ ', a1)
        commandline = paste0('use_labels(data2work, glm(formula(', tmp_formula_str, '), family = binomial(link = "probit"), data = data2work ) )' )
        tmp_model = eval(parse(text = commandline))
        # cov1         = vcovHC(tmp_model, type = "HC1") # Never do robust errors for non linear model. Sometimes they are not even consistent. 
        # robust_se[[i]] = sqrt(diag(cov1))
        reg_list[[i]] = tmp_model
        x_vars = c(x_vars, a2)
        i = i + 1
    }
}

res = remove_backticks(
    stargazer(
        reg_list,
        title="Results",
        align=TRUE,
        omit.stat = c("f"),
        dep.var.caption  =  paste0("Dependent variable: ", eval(parse(text = paste0('var_lab(data2work$',dep_list[1],')'))) ),
        dep.var.labels.include = FALSE,
        df = FALSE,
        digits = 5,
        # notes        = TeX("Sometimes you just have to start over."),
        # notes.append = FALSE,
        se = robust_se))
# res[grepl("Note",res)] = note.latex
res[grepl("\\label\\{\\}",res)] = "%"
res[grepl("\\caption\\{",res)] = "%"
res[grepl("Table created by stargazer v.5.2.2 by Marek Hlavac",res)] = ""
res[grepl("\\% Date and time:",res)] = ""
res[grepl("\\% Requires LaTeX packages: dcolumn",res)] = ""
regmatches(res, gregexpr("(-?)(\\d+(,\\d+)?(\\.\\d+)?)", res)) =  lapply(regmatches(res, gregexpr("(-?)(\\d+(,\\d+)?(\\.\\d+)?)", res) ) , mysignif )
# 
res[grepl("begin\\{table\\}",res)] = begintable.latex
res[grepl("end\\{table\\}",res)] = endtable.latex
cat("\014")  

write_clip(res)
writeLines(res)
cat("\014")

# End



tmp_model = glm( critical_flag ~ tdiff50_recv1 + Repo_balances1 + coupon_issuances + QuarterEnd, family = binomial(link = "probit"), data = data2work )

tem1 = sd(data2work$tdiff50_recv1, na.rm=TRUE)
tem2 = sd(data2work$Repo_balances1, na.rm=TRUE)
# tem3 = sd(data2work$coupon_issuances, na.rm=TRUE)
tem3 = mean(data2work$coupon_issuances[data2work$coupon_issuances > 0.01], method = 'mfv', na.rm=TRUE)

# dataworktem = data2work %>% select(coupon_issuances) %>% filter(coupon_issuances > 0.01) %>% mutate(meancoupon = mean(coupon_issuances) )

tem4 = mean(data2work$tdiff50_recv1, na.rm=TRUE)
tem5 = mean(data2work$Repo_balances1, na.rm=TRUE)
tem6 = mean(data2work$coupon_issuances, na.rm=TRUE)

data2work_new = tibble(
    tdiff50_recv1	    =	c(	tem4	,	tem4+tem1	,	 tem4	    ,  tem4+tem1, tem4+tem1  ,   	 tem4	   ,  	 tem4+tem1 ,	tem4	,	tem4+tem1	,	 tem4	     ,	 tem4	),
    Repo_balances1	    =	c(	tem5	,	tem5	    ,	 tem5-tem2	,  tem5-tem2, tem5-tem2  ,   	 tem5	   ,  	 tem5-tem2 ,	tem5	,	tem5	    ,	 tem5-tem2	 ,	 tem5	),
    coupon_issuances	=	c(	tem6	,	tem6	    ,	 tem6	    ,  tem6	    , tem6+tem3	 ,   	 tem6+tem3 ,  	 tem6+tem3 ,	tem6	,	tem6	    ,	 tem6	     ,	 tem6+tem3	),
    QuarterEnd 	        =	c(	0	    ,	0	        ,	0	        ,   0	    ,  0	     ,   	0	       ,  	0	       ,	1	    ,	1	        ,	1	         ,	1	)
    # tdiff50_recv1	=	c(	tem4	,	tem4+tem1	),
    # Repo_balances1	=	c(	tem5	,	 tem5+tem2	),
    # coupon_issuances	=	c(	 tem6	,	 tem6+tem3	),
    # QuarterEnd 	=	c(	0	,	1	)		
)

result = predict(tmp_model, type = "response", newdata = data2work_new, se.fit = TRUE)

data2work_new$prediction = as.numeric(result$fit) 
data2work_new$se = as.numeric(result$se.fit) 

data2work_new %<>% select(prediction, se, tdiff50_recv1, Repo_balances1, coupon_issuances,  QuarterEnd ) %>%
    mutate(Repo_balances1 = Repo_balances1*1000, coupon_issuances = coupon_issuances*1000 )

res = stargazer(data2work_new, summary=FALSE, rownames=FALSE, digits = 4)

# Number inputs for Figure 5
regmatches(res, gregexpr("(-?)(\\d+(\\.\\d+)?)", res)) =  lapply(regmatches(res, gregexpr("(-?)(\\d+(\\.\\d+)?)", res) ) , mysignif )

write_clip(res)
writeLines(res)

tem1 = sd(data2work$tdiff50_recv1, na.rm=TRUE)
tem2 = sd(data2work$Repo_balances1, na.rm=TRUE)
tem3 = sd(data2work$coupon_issuances, na.rm=TRUE)

tem4 = mean(data2work$tdiff50_recv1, na.rm=TRUE)
tem5 = mean(data2work$Repo_balances1, na.rm=TRUE)
tem6 = mean(data2work$coupon_issuances, na.rm=TRUE)
c(tem1, 
  tem2, 
  tem3, 
  tem4, 
  tem5, 
  tem6)

# End




# IV regression---------------------------
# For all IV results, you need to start from here and do not jump around. 

rm( list = setdiff( ls(), c(keepVariables, "") ))
cat("\014")  
data2work = data2work1
# data2work = data2work1 %>% filter(Date <= ymd('2020-10-30'))


# ------------Stage 1 regression---------------------------------------------------------------------------------------------------------------------------------------------
# Table 13
model = use_labels(data2work, 
                   lm(Repo_balances1 ~ lag(p25Diff,n = 7) + lag(p50Diff,n = 7) + lag(p75Diff,n = 7) + 
                          Tbills_outstanding + total_redemptions_all + 
                          total_issuance_bills + coupon_issuances + QuarterEnd + DlrDepUnins + Corporation_Income_Taxes_trillion, 
                      na.action=na.exclude,
                      data = data2work ) )


model.1 = use_labels(data2work, 
                     lm(
                         tdiff50_recv1 ~ lag(Repo_balances1),
                         # Repo_balances1 ~ lag(tdiff50_recv1),
                         na.action=na.exclude,
                         data = data2work ) )

model.2 = use_labels(data2work, 
                     lm(
                         # tdiff50_recv1 ~ lag(Repo_balances1),
                         Repo_balances1 ~ lag(tdiff50_recv1),
                         na.action=na.exclude,
                         data = data2work ) )

tem = data2work
tem$Repo_balance_hat = fitted(model)
tem$Repo_balance_residual = residuals(model.2)
tem$tdiff50_recv1_residual = residuals(model.1)

tem %<>% select(Date, Repo_balance_hat, Repo_balance_residual, tdiff50_recv1_residual)

data2work7 = left_join(data2work, tem, by = "Date")
data2work7 =
    apply_labels(
        data2work7,
        Repo_balance_hat = "predicted dealer opening balances", 
        Repo_balance_residual = "residual of dealer opening balances on lag median receive time from all banks", 
        tdiff50_recv1_residual = "residual of median receive time on opening balances"
    )
rm( tem )

res = remove_backticks( 
    stargazer(
        model,
        title="Results", 
        align=TRUE, 
        # omit.stat = c("f"), 
        dep.var.caption  =  paste0("Dependent variable: ", eval(parse(text = paste0('var_lab(data2work$', 'Repo_balances1',')'))) ),
        dep.var.labels.include = FALSE,
        df = FALSE,
        digits = 5
    ))
res[grepl("\\label\\{\\}",res)] = "%"
res[grepl("\\caption\\{",res)] = "%"

res[grepl("Table created by stargazer v.",res)] = ""
res[grepl("\\% Date and time:",res)] = ""
res[grepl("\\% Requires LaTeX packages: dcolumn",res)] = ""

regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res)) =  lapply(regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res) ) , mysignif )
res[grepl("begin\\{table\\}",res)] = begintable.latex
res[grepl("end\\{table\\}",res)] = endtable.latex
cat("\014")  
write_clip(res)
writeLines(res)
cat("\014")  


# ---------------------------------Stage 2 regression-------------------------------------------------------------------------------------------------------------------------------------------
# Table 12

library("ivreg")
# library("AER")


regressor_list = c(
    'Repo_balances1 + tdiff50_recv1 + QuarterEnd + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins + Corporation_Income_Taxes_trillion', 
    # 
    'Repo_balances1 + QuarterEnd + Tbills_outstanding  + total_redemptions_all +
     total_issuance_bills + coupon_issuances + DlrDepUnins + Corporation_Income_Taxes_trillion | lag(p25Diff,n = 7) +
     lag(p50Diff,n = 7) + lag(p75Diff,n = 7) + Tbills_outstanding + total_redemptions_all +
                          total_issuance_bills + coupon_issuances + QuarterEnd + DlrDepUnins + Corporation_Income_Taxes_trillion',
    # 
    # #
    'Repo_balances1 + tdiff50_recv1_residual + QuarterEnd + Tbills_outstanding  + total_redemptions_all +
     total_issuance_bills + coupon_issuances + DlrDepUnins + Corporation_Income_Taxes_trillion | lag(p25Diff,n = 7) +
     lag(p50Diff,n = 7) + lag(p75Diff,n = 7) + tdiff50_recv1_residual + Tbills_outstanding + total_redemptions_all +
                          total_issuance_bills + coupon_issuances + QuarterEnd + DlrDepUnins + Corporation_Income_Taxes_trillion'
    # #
    )



dep_list = c('SOFR_IOER')

reg_list = vector(mode = 'list', length = length(regressor_list) * length(dep_list) )
robust_se = vector(mode = 'list', length = length(regressor_list) * length(dep_list)  )
x_vars = c()
i = 1

for (a1 in regressor_list) {
    for (a2 in dep_list) {
        tmp_formula_str = paste0(a2, ' ~ ', a1)
        commandline = paste0('use_labels(data2work7, ivreg(formula(', tmp_formula_str, '), data = data2work7 ) )' )
        tmp_model = eval(parse(text = commandline))
        cov1         = vcovHC(tmp_model, type = "HC1")
        robust_se[[i]] = sqrt(diag(cov1))
        reg_list[[i]] = tmp_model
        x_vars = c(x_vars, a2)
        i = i + 1
    }
}




note.latex = paste0(
    "\\multicolumn{",
    length(reg_list)+1,
    "}{l}{{Note:} Standard errors are adjusted for heteroskedasticity. $^{*}p<0.1$; $^{**}p<0.05$; $^{***}p<0.01$.} \\\\ \\multicolumn{",
    length(reg_list)+1,
    "}{l}{ \\qquad \\quad Constant included for each specification.}")

res = remove_backticks(
    stargazer(
        reg_list,
        title="Results",
        align=TRUE,
        omit.stat = c("f"),
        dep.var.caption  =  paste0("Dependent variable: ", eval(parse(text = paste0('var_lab(data2work$',dep_list[1],')'))) ),
        dep.var.labels.include = FALSE,
        df = FALSE,
        digits = 5,
        # notes        = TeX("Sometimes you just have to start over."),
        # notes.append = FALSE,
        se = robust_se
        ))
res[grepl("Note",res)] = note.latex
res[grepl("\\label\\{\\}",res)] = "%"
res[grepl("\\caption\\{",res)] = "%"
res[grepl("Table created by stargazer v.",res)] = ""
res[grepl("\\% Date and time:",res)] = ""
res[grepl("\\% Requires LaTeX packages: dcolumn",res)] = ""
regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res)) =  lapply(regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res) ) , mysignif )
res[grepl("begin\\{table\\}",res)] = begintable.latex
res[grepl("end\\{table\\}",res)] = endtable.latex
cat("\014")  

write_clip(res)
writeLines(res)
# End



# # ---------IV Correcting STD
rm( list = setdiff( ls(), c(keepVariables, "") ))
cat("\014")  
data2work = data2work1

data2work = data2work1 %>% arrange(Date)  # Use first definition
# ------------Stage 1 regression

model = use_labels(data2work, 
                   lm(Repo_balances1 ~ lag(p25Diff,n = 7) + lag(p50Diff,n = 7) + lag(p75Diff,n = 7) + 
                          Tbills_outstanding + total_redemptions_all + 
                          total_issuance_bills + coupon_issuances + QuarterEnd + DlrDepUnins + Corporation_Income_Taxes_trillion, 
                      na.action=na.exclude,
                      data = data2work ) )

model.0 = lm(Repo_balances1 ~ lag(p25Diff,n = 7) + lag(p50Diff,n = 7) + lag(p75Diff,n = 7) + 
                 Tbills_outstanding + total_redemptions_all + 
                 total_issuance_bills + coupon_issuances + QuarterEnd + DlrDepUnins + Corporation_Income_Taxes_trillion, 
             na.action=na.exclude,
             data = data2work )

model.1 = use_labels(data2work, 
                     lm(
                         # tdiff50_recv1 ~ lag(Repo_balances1),
                         tdiff50_recv1 ~ Repo_balances1,
                         na.action=na.exclude,
                         data = data2work ) )



tem = data2work
tem$Repo_balance_hat = fitted(model)
tem$tdiff50_recv1_residual = residuals(model.1)

tem %<>% select(Date, Repo_balance_hat, tdiff50_recv1_residual)

data2work7 = left_join(data2work, tem, by = "Date")
data2work7 =
    apply_labels(
        data2work7,
        Repo_balance_hat = "predicted dealer opening balances", 
        # Repo_balance_residual = "residual of dealer opening balances on lag median receive time from all banks", 
        tdiff50_recv1_residual = "residual of median receive time on opening balances"
    )
rm( tem )

data2work_tem = data2work %>% 
    mutate(
        # error = Repo_balances1 - model$coefficients[1] - lag(tdiff25_recvF1)* model$coefficients[2] - lag(tdiff50_recvF1)*model$coefficients[3] -
        #     lag(tdiff75_recvF1)*model$coefficients[4] - Corporation_Income_Taxes_trillion*model$coefficients[5],
        error =Repo_balances1 - ( model$coefficients[1] +
            lag(p25Diff,n = 7)*model$coefficients[2] + lag(p50Diff,n = 7)*model$coefficients[3] + lag(p75Diff,n = 7)*model$coefficients[4] + 
            Tbills_outstanding*model$coefficients[5] + total_redemptions_all*model$coefficients[6] + 
            total_issuance_bills*model$coefficients[7] + coupon_issuances*model$coefficients[8] + QuarterEnd*model$coefficients[9] + DlrDepUnins*model$coefficients[10] + Corporation_Income_Taxes_trillion*model$coefficients[11]),
        error2 = tdiff50_recv1 - Repo_balances1 * model.1$coefficients[2] - model.1$coefficients[1] ) %>% 
    mutate(
        lag_tdiff25_recvF1   = lag(tdiff25_recvF1),  
        lag_tdiff50_recvF1   = lag(tdiff50_recvF1),    
        lag_tdiff75_recvF1   = lag(tdiff75_recvF1),
        lag_p25Diff = lag(p25Diff,n = 7),  
        lag_p50Diff = lag(p50Diff,n = 7), 
        lag_p75Diff = lag(p75Diff,n = 7) 
        ) %>% 
    mutate(covarerror = error*error2) %>% filter(!is.na(covarerror) )

tem = data2work_tem %>% select(Date, lag_p25Diff, lag_p50Diff, lag_p75Diff)

data2work7 = left_join(data2work7, tem, by = "Date")
rm( tem )




Xtem = cbind(
    rep(1, 2027),
    data2work_tem$lag_p25Diff,
    data2work_tem$lag_p50Diff,
    data2work_tem$lag_p75Diff,
    data2work_tem$Tbills_outstanding,
    data2work_tem$total_redemptions_all,
    data2work_tem$total_issuance_bills,
    data2work_tem$coupon_issuances,
    data2work_tem$QuarterEnd,
    data2work_tem$DlrDepUnins,
    data2work_tem$Corporation_Income_Taxes_trillion
)
Ytem = cbind(
    rep(1, 2027),
    data2work_tem$Repo_balances1
)
covtemvec = solve(t(Xtem) %*% Xtem ) %*% t(Xtem) %*% diag(data2work_tem$covarerror) %*% Ytem %*% solve(t(Ytem) %*% Ytem ) 
# model.1$residuals
vcov(model)
vcov(model.1)
# Testing assuming no correlation
# Matrix_tem1 = cbind( vcov(model), c(0,0,0,0,0),  c(0,0,0,0,0))
# Matrix_tem2 = cbind( c(0,0),  c(0,0), c(0,0), c(0,0), c(0,0), vcov(model.1))
# Matrix_tem3 = rbind(Matrix_tem1, Matrix_tem2)
Matrix_tem1 = cbind( vcov(model), covtemvec )
Matrix_tem2 = cbind( t(covtemvec), vcov(model.1))
Matrix_tem = rbind(Matrix_tem1, Matrix_tem2)
# mean(residuals(model.1) * residuals(model), na.rm = TRUE) = 1.90 assume it is zero to simplify my calculation. 
# sd(residuals(model.1), na.rm = TRUE)
# sd(residuals(model), na.rm = TRUE)
# mean(residuals(model.1), na.rm = TRUE)
# mean(residuals(model), na.rm = TRUE)
sd(residuals(model.1), na.rm = TRUE) * sd(residuals(model), na.rm = TRUE)
1.90/(sd(residuals(model.1), na.rm = TRUE) * sd(residuals(model), na.rm = TRUE))

# ---------------------------------Stage 2 regression, correcting std for IV
library("ivreg")

# We need to adjust the standard deviation since the variable 
# residual of median receive time on opening balances
# is generated from model.2. 
# Here the adjustment comes application of Kevin M. Murphy, Robert H. Topel (1985)

regressor_list = c(
    'Repo_balances1 + tdiff50_recv1_residual + QuarterEnd + Tbills_outstanding  + total_redemptions_all +
     total_issuance_bills + coupon_issuances + DlrDepUnins + Corporation_Income_Taxes_trillion | lag(p25Diff,n = 7) +
     lag(p50Diff,n = 7) + lag(p75Diff,n = 7) + tdiff50_recv1_residual + Tbills_outstanding + total_redemptions_all +
                          total_issuance_bills + coupon_issuances + QuarterEnd + DlrDepUnins + Corporation_Income_Taxes_trillion'
)
dep_list = c('SOFR_IOER')
# dep_list = c('GCF_IOER')
# dep_list = c('TGCR_IOER')

reg_list = vector(mode = 'list', length = length(regressor_list) * length(dep_list) )
robust_se = vector(mode = 'list', length = length(regressor_list) * length(dep_list)  )
x_vars = c()
i = 1

for (a1 in regressor_list) {
    for (a2 in dep_list) {
        tmp_formula_str = paste0(a2, ' ~ ', a1)
        commandline = paste0('use_labels(data2work7, ivreg(formula(', tmp_formula_str, '), data = data2work7 ) )' )
        tmp_model = eval(parse(text = commandline))
        cov1         = vcovHC(tmp_model, type = "HC1")
        robust_se[[i]] = sqrt(diag(cov1))
        reg_list[[i]] = tmp_model
        x_vars = c(x_vars, a2)
        i = i + 1
    }
}


note.latex = paste0(
    "\\multicolumn{",
    length(reg_list)+1,
    "}{l}{{Note:} Standard errors are adjusted for heteroskedasticity. $^{*}p<0.1$; $^{**}p<0.05$; $^{***}p<0.01$.} \\\\ \\multicolumn{",
    length(reg_list)+1,
    "}{l}{ \\qquad \\quad Constant included for each specification.}")

res = remove_backticks(
    stargazer(
        reg_list,
        title="Results",
        align=TRUE,
        omit.stat = c("f"),
        dep.var.caption  =  paste0("Dependent variable: ", eval(parse(text = paste0('var_lab(data2work$',dep_list[1],')'))) ),
        dep.var.labels.include = FALSE,
        df = FALSE,
        digits = 5,
        # notes        = TeX("Sometimes you just have to start over."),
        # notes.append = FALSE,
        se = robust_se))
res[grepl("Note",res)] = note.latex
res[grepl("\\label\\{\\}",res)] = "%"
res[grepl("\\caption\\{",res)] = "%"
res[grepl("Table created by stargazer v.5.2.2 by Marek Hlavac",res)] = ""
res[grepl("\\% Date and time:",res)] = ""
res[grepl("\\% Requires LaTeX packages: dcolumn",res)] = ""
regmatches(res, gregexpr("(-?)(\\d+(\\.\\d+)?)", res)) =  lapply(regmatches(res, gregexpr("(-?)(\\d+(\\.\\d+)?)", res) ) , mysignif )
# 
# 
res[grepl("begin\\{table\\}",res)] = begintable.latex
res[grepl("end\\{table\\}",res)] = endtable.latex
cat("\014")  

write_clip(res)
writeLines(res)



model2work = reg_list[[1]]




Z = cbind( rep(1,2061),  data2work7$QuarterEnd,  data2work7$Tbills_outstanding, data2work7$total_redemptions_all,  
           data2work7$total_issuance_bills, data2work7$coupon_issuances, data2work7$DlrDepUnins,  data2work7$Corporation_Income_Taxes_trillion, data2work7$Repo_balance_hat,  data2work7$tdiff50_recv1_residual)



Fstar_tem1 = model2work$coefficients[[2]] * cbind( rep(1,2061),  
                                                   data2work7$lag_p25Diff,  
                                                   data2work7$lag_p50Diff,   
                                                   data2work7$lag_p75Diff,  
                                                   data2work7$Tbills_outstanding,
                                                   data2work7$total_redemptions_all,
                                                   data2work7$total_issuance_bills,
                                                   data2work7$coupon_issuances,
                                                   data2work7$QuarterEnd,
                                                   data2work7$DlrDepUnins,
                                                   data2work7$Corporation_Income_Taxes_trillion) 

Fstar_tem2 = - model2work$coefficients[[3]] * cbind( rep(1,2061),  data2work7$Repo_balances1) 
Fstar = cbind(Fstar_tem1, Fstar_tem2)
Vec_tem = cbind(Z, Fstar)
Vec_tem <- Vec_tem[complete.cases(Vec_tem), ]  

Fstar <- Vec_tem[,11:23]  
Z <- Vec_tem[,1:10]  


A = t(Z)%*%Z
B = t(Z)%*%Fstar

result_tem = solve(A) %*% B %*%  Matrix_tem %*% t(B) %*%  (solve(A)) 

result_vec_tem = diag(result_tem)
result_vec_tem  + c(2.53,
                    3.18, 
                    0.863,
                    7.31,
                    6.88,
                    6.39,
                    3.53,
                    33.9,
                    8.91,
                    0.0272)
# adding results from the previous results. 
# 3.1960615 (this constant)   3.5266746   0.8949359  21.8276714  20.5183734  16.7988902   5.1673044 198.7552320  13.7902650   0.0272006 (this last two is the first two)

23.3/13.79


                    
2*(1 - pnorm(46.6, 0, 8.15  ) )












# ----------------------Quantile Regression -----------March 2023--------------------------
library(quantreg)
rm( list = setdiff( ls(), c(keepVariables, "") ))
cat("\014")  
data2work = data2work1

# ----------------------Main specification, no IV; 99% level--------------------
# Table 8 
library(quantreg)
data2work = data2work1

rm( list = setdiff( ls(), c(keepVariables, "") ))
cat("\014")

regressor_list = c(
    'Repo_balances1', 
    'Repo_balances1 + QuarterEnd',
    'tdiff50_recv1  + QuarterEnd',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + DlrDepUnins',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins + Corporation_Income_Taxes_trillion'
)


dep_list = c('SOFR_IOER')

reg_list = vector(mode = 'list', length = length(regressor_list) * length(dep_list) )
robust_se = vector(mode = 'list', length = length(regressor_list) * length(dep_list) )
Rsquare_list = vector(mode = 'character', length = length(regressor_list) * length(dep_list) )

tmp_model0 = rq(SOFR_IOER ~ 1, data = data2work,  tau = 0.99)

x_vars = c()
i = 1
for (a1 in regressor_list) {
    for (a2 in dep_list) {
        tmp_formula_str = paste0(a2, ' ~ ', a1)
        commandline = paste0('use_labels(data2work, rq(formula(', tmp_formula_str, '), data = data2work,  tau = 0.99) )' )
        tmp_model = eval(parse(text = commandline))
        reg_list[[i]] = tmp_model
        # cov1         = vcovHC(tmp_model, type = "HC1")
        # robust_se[[i]] = sqrt(diag(cov1))
        Rsquare_list[[i]] = 1 - tmp_model$rho/tmp_model0$rho
        x_vars = c(x_vars, a2)
        i = i + 1
    }
}


note.latex = paste0(
    "\\multicolumn{",
    length(reg_list)+1,
    "}{l}{{Note:} Constant included for each specification.}")


res = remove_backticks(
    stargazer(
        reg_list,
        title="Results",
        align=TRUE,
        omit.stat = c("f"),
        dep.var.caption  =  paste0("Dependent variable: ", eval(parse(text = paste0('var_lab(data2work$',dep_list[1],')'))) ),
        dep.var.labels.include = FALSE,
        df = FALSE,
        digits = 5,
        # notes        = TeX("Sometimes you just have to start over."),
        # notes.append = FALSE,
        # se = robust_se
        add.lines=list( c("pseudo-$R^{1}$", Rsquare_list) ), 
        table.layout = "=lc#-t-sa-n"
        )
    )
res[grepl("Note",res)] = note.latex
res[grepl("\\label\\{\\}",res)] = "%"
res[grepl("\\caption\\{",res)] = "%"
res[grepl("Table created by stargazer v.5.2.3 by Marek Hlavac",res)] = ""
res[grepl("\\% Date and time:",res)] = ""
res[grepl("\\% Requires LaTeX packages: dcolumn",res)] = ""
regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res)) =  lapply(regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res) ) , mysignif )
res[grepl("begin\\{table\\}",res)] = begintable.latex
res[grepl("end\\{table\\}",res)] = endtable.latex
cat("\014")  
write_clip(res)
writeLines(res)
# End



# ----------------------Table 2 99% level--------------------
library(quantreg)
data2work = data2work1

# Table 2 Column 1

rm( list = setdiff( ls(), c(keepVariables, "") ))
cat("\014")

regressor_list = c(
    'Repo_balances1 + QuarterEnd + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins + Corporation_Income_Taxes_trillion',
    'Repo_balances1 + QuarterEnd + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins + Corporation_Income_Taxes_trillion'
)


# dep_list = c('SOFR_IOER', 'GCF_IOER')
dep_list = c('SOFR_IOER')

reg_list = vector(mode = 'list', length = length(regressor_list) * length(dep_list) )
robust_se = vector(mode = 'list', length = length(regressor_list) * length(dep_list) )
Rsquare_list = vector(mode = 'character', length = length(regressor_list) * length(dep_list) )

tmp_model0 = rq(SOFR_IOER ~ 1, data = data2work,  tau = 0.99)

x_vars = c()
i = 1
for (a2 in dep_list) {
for (a1 in regressor_list) {
        tmp_formula_str = paste0(a2, ' ~ ', a1)
        commandline = paste0('use_labels(data2work, rq(formula(', tmp_formula_str, '), data = data2work,  tau = 0.99) )' )
        tmp_model = eval(parse(text = commandline))
        reg_list[[i]] = tmp_model
        # cov1         = vcovHC(tmp_model, type = "HC1")
        # robust_se[[i]] = sqrt(diag(cov1))
        Rsquare_list[[i]] = 1 - tmp_model$rho/tmp_model0$rho
        x_vars = c(x_vars, a2)
        i = i + 1
    }
}

note.latex = paste0(
    "\\multicolumn{",
    length(reg_list)+1,
    "}{l}{{Note:} Constant included for each specification.}")


res = remove_backticks(
    stargazer(
        reg_list,
        title="Results",
        align=TRUE,
        omit.stat = c("f"),
        dep.var.caption  =  paste0("Dependent variable: ", eval(parse(text = paste0('var_lab(data2work$',dep_list[1],')'))) ),
        dep.var.labels.include = FALSE,
        df = FALSE,
        digits = 5,
        # notes        = TeX("Sometimes you just have to start over."),
        # notes.append = FALSE,
        # se = robust_se
        add.lines=list( c("pseudo-$R^{1}$", Rsquare_list) ),
        table.layout = "=lc#-t-sa-n"
    )
)
res[grepl("Note",res)] = note.latex
res[grepl("\\label\\{\\}",res)] = "%"
res[grepl("\\caption\\{",res)] = "%"
res[grepl("Table created by stargazer v.5.2.3 by Marek Hlavac",res)] = ""
res[grepl("\\% Date and time:",res)] = ""
res[grepl("\\% Requires LaTeX packages: dcolumn",res)] = ""
regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res)) =  lapply(regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res) ) , mysignif )
res[grepl("begin\\{table\\}",res)] = begintable.latex
res[grepl("end\\{table\\}",res)] = endtable.latex
cat("\014")
write_clip(res)
writeLines(res)


# To calculate the fake R1 for the IV quantile regression. (Estimates are generated from the stata file. )
# Caution: the R1 for IV quantile is not well studied, and hence not well defined in the literature. Here we apply the algorithm of R1 as in standard 
# quantile. However just like R-square in 2 stage linear regression, the R1 here is hard to interpret so readers who want to replicate should understand
# the danger of using such R1. In the paper we do not discuss the meaning of R1; We compute the number just to complete the table. 

weights <- tmp_model$coef
weights2 = c(6.34	,
             -74.1	,
             31.5	,
             3.22	,
             -13.7	,
             23.2	,
             38.8	,
             24.1	,
             15.2	)
for (index in c(1:length(weights2))) {
    weights[index] = weights2[index]
}

datatem <- model.frame(formula = formula(tmp_model), data = model.frame(tmp_model),na.action=NULL)
datatem_X = datatem %>% select(-`SOFR - IOR`)
datatem_Y = tmp_model$y
# mm <- model.matrix(formula(tmp_model, fixed.only = TRUE), data = datatem)
# offset <- as.vector(mm %*% weights)
# resid = datatem$`SOFR - IOR` - offset
# rho_IVQuantile = 0.99 * pmax(resid, 0) - (1- 0.99) * pmin(resid, 0)
# V_IVQuantile = sum(rho_IVQuantile)
# rho_IVQuantile2 = 0.99 * pmax(tmp_model0$residuals, 0) - (1- 0.99) * pmin(tmp_model0$residuals, 0)
# V_IVQuantile2 = sum(rho_IVQuantile2)

tmp_model_IV = tmp_model
tmp_model_IV$coefficients = weights   # Be very careful here. This only only works in this context. 
# tmp_model$coefficients # double checking
residual_IV = datatem_Y - predict(tmp_model_IV, newdata = datatem_X, type = c("response"))
tmp_model_IV$rho = sum(residual_IV*(0.99-(residual_IV<0)))
R1_IV = 1 - tmp_model_IV$rho/tmp_model0$rho
R1_IV
# End

# ----------------------Quantile Regression -----------March 2023--------------------------
# ----------------------Main specification, no IV; 99% level No Time Receive--------------------
library(quantreg)
data2work = data2work1

rm( list = setdiff( ls(), c(keepVariables, "") ))
cat("\014")
regressor_list = c(
    'Repo_balances1',
    'Repo_balances1 + QuarterEnd',
    'tdiff50_recv1  + QuarterEnd',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + DlrDepUnins',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins',
    # 'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + NetPosition_try + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins + Corporation_Income_Taxes_trillion'
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins + Corporation_Income_Taxes_trillion'
)

dep_list = c('SOFR_IOER')
reg_list = vector(mode = 'list', length = length(regressor_list) * length(dep_list) )
robust_se = vector(mode = 'list', length = length(regressor_list) * length(dep_list) )
Rsquare_list = vector(mode = 'character', length = length(regressor_list) * length(dep_list) )

tmp_model0 = rq(SOFR_IOER ~ 1, data = data2work,  tau = 0.99)

x_vars = c()
i = 1
for (a1 in regressor_list) {
    for (a2 in dep_list) {
        tmp_formula_str = paste0(a2, ' ~ ', a1)
        commandline = paste0('use_labels(data2work, rq(formula(', tmp_formula_str, '), data = data2work,  tau = 0.99) )' )
        tmp_model = eval(parse(text = commandline))
        reg_list[[i]] = tmp_model
        # cov1         = vcovHC(tmp_model, type = "HC1")
        # robust_se[[i]] = sqrt(diag(cov1))
        Rsquare_list[[i]] = 1 - tmp_model$rho/tmp_model0$rho
        x_vars = c(x_vars, a2)
        i = i + 1
    }
}

note.latex = paste0(
    "\\multicolumn{",
    length(reg_list)+1,
    "}{l}{{Note:} Constant included for each specification.}")


res = remove_backticks(
    stargazer(
        reg_list,
        title="Results",
        align=TRUE,
        omit.stat = c("f"),
        dep.var.caption  =  paste0("Dependent variable: ", eval(parse(text = paste0('var_lab(data2work$',dep_list[1],')'))) ),
        dep.var.labels.include = FALSE,
        df = FALSE,
        digits = 5,
        # notes        = TeX("Sometimes you just have to start over."),
        # notes.append = FALSE,
        # se = robust_se
        add.lines=list( c("pseudo-$R^{1}$", Rsquare_list) ), 
        table.layout = "=lc#-t-sa-n"
    )
)
res[grepl("Note",res)] = note.latex
res[grepl("\\label\\{\\}",res)] = "%"
res[grepl("\\caption\\{",res)] = "%"
res[grepl("Table created by stargazer v.5.2.3 by Marek Hlavac",res)] = ""
res[grepl("\\% Date and time:",res)] = ""
res[grepl("\\% Requires LaTeX packages: dcolumn",res)] = ""
regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res)) =  lapply(regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res) ) , mysignif )
res[grepl("begin\\{table\\}",res)] = begintable.latex
res[grepl("end\\{table\\}",res)] = endtable.latex
cat("\014")  
write_clip(res)
writeLines(res)
# End


# ----------------------Main specification, no IV; 95% level--------------------
# Table 7
rm( list = setdiff( ls(), c(keepVariables, "") ))
cat("\014")
data2work = data2work1

rm( list = setdiff( ls(), c(keepVariables, "") ))
cat("\014")
regressor_list = c(
    'Repo_balances1',
    'Repo_balances1 + QuarterEnd',
    'tdiff50_recv1  + QuarterEnd',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + DlrDepUnins',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins',
    'Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins + Corporation_Income_Taxes_trillion'
)

dep_list = c('SOFR_IOER')
# dep_list = c('GCF_IOER')
# dep_list = c('TGCR_IOER')

reg_list = vector(mode = 'list', length = length(regressor_list) * length(dep_list) )
robust_se = vector(mode = 'list', length = length(regressor_list) * length(dep_list) )
Rsquare_list = vector(mode = 'character', length = length(regressor_list) * length(dep_list) )

tmp_model0 = rq(SOFR_IOER ~ 1, data = data2work,  tau = 0.95)

# rho <- function(u,tau=.5)u*(tau - (u < 0))
# # R1 <- 1 - fit1$rho/fit0$rho
# sum(rho(tmp_model0$resid, tmp_model0$tau))

x_vars = c()
i = 1
for (a1 in regressor_list) {
    for (a2 in dep_list) {
        tmp_formula_str = paste0(a2, ' ~ ', a1)
        commandline = paste0('use_labels(data2work, rq(formula(', tmp_formula_str, '), data = data2work,  tau = 0.95) )' )
        tmp_model = eval(parse(text = commandline))
        reg_list[[i]] = tmp_model
        # cov1         = vcovHC(tmp_model, type = "HC1")
        # robust_se[[i]] = sqrt(diag(cov1))
        Rsquare_list[[i]] = 1 - tmp_model$rho/tmp_model0$rho
        x_vars = c(x_vars, a2)
        i = i + 1
    }
}


note.latex = paste0(
    "\\multicolumn{",
    length(reg_list)+1,
    "}{l}{{Note:} Constant included for each specification.}")

res = remove_backticks(
    stargazer(
        reg_list,
        title="Results",
        align=TRUE,
        omit.stat = c("f"),
        dep.var.caption  =  paste0("Dependent variable: ", eval(parse(text = paste0('var_lab(data2work$',dep_list[1],')'))) ),
        dep.var.labels.include = FALSE,
        df = FALSE,
        digits = 5,
        # notes        = TeX("Sometimes you just have to start over."),
        # notes.append = FALSE,
        # se = robust_se
        add.lines=list( c("pseudo-$R^{1}$", Rsquare_list) ), 
        table.layout = "=lc#-t-sa-n"
    )
)
res[grepl("Note",res)] = note.latex
res[grepl("\\label\\{\\}",res)] = "%"
res[grepl("\\caption\\{",res)] = "%"
res[grepl("Table created by stargazer v.5.2.3 by Marek Hlavac",res)] = ""
res[grepl("\\% Date and time:",res)] = ""
res[grepl("\\% Requires LaTeX packages: dcolumn",res)] = ""
regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res)) =  lapply(regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res) ) , mysignif )
res[grepl("begin\\{table\\}",res)] = begintable.latex
res[grepl("end\\{table\\}",res)] = endtable.latex
cat("\014")  
write_clip(res)
writeLines(res)
# End

tmp_model0 = rq(SOFR_IOER ~ Repo_balances1 + QuarterEnd + tdiff50_recv1, data = data2work,  tau = 0.95)
tmp_model1 = rq(SOFR_IOER ~ Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + NetPosition_try + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins, data = data2work,  tau = 0.95)
1 - tmp_model$rho/tmp_model0$rho



# -----------------Ploting------------------------------------
library(quantreg)

data2work = data2work1
rm( list = setdiff( ls(), c(keepVariables, "") ))
cat("\014")

tau_vector = seq(from = 0.51, to = 0.99, by = 0.01)

coeff_vector_dealer = 0 * tau_vector
coeff_vector_time = 0 * tau_vector
se_vector_dealer= 0 * tau_vector
se_vector_time = 0 * tau_vector

coeff_vector_dealer2 = 0 * tau_vector
coeff_vector_time2 = 0 * tau_vector
se_vector_dealer2= 0 * tau_vector
se_vector_time2 = 0 * tau_vector

reg_list = vector(mode = 'list', length = length(tau_vector) )
robust_se = vector(mode = 'list', length = length(tau_vector) )
Rsquare_list = 0 * tau_vector

reg_list2 = vector(mode = 'list', length = length(tau_vector) )
robust_se2 = vector(mode = 'list', length = length(tau_vector) )
Rsquare_list2 = 0 * tau_vector



i = 0
for (tau_tem in tau_vector){
    print(tau_tem)
    i = i + 1
    tmp_model0 = rq(SOFR_IOER ~ 1, data = data2work,  tau = tau_tem)
    tmp_formula_str = paste0('SOFR_IOER ~ Repo_balances1 + QuarterEnd + tdiff50_recv1 + DlrDepUnins')
    commandline = paste0('use_labels(data2work, rq(formula(', tmp_formula_str, '), data = data2work,  tau = ', tau_tem, ') )' )
    tmp_model = eval(parse(text = commandline))
    coeff_vector_dealer[i] = tmp_model$coefficients[2]
    coeff_vector_time[i] = tmp_model$coefficients[4]
    reg_list[[i]] = tmp_model
    Rsquare_list[i] = 1 - tmp_model$rho/tmp_model0$rho
    
    
    # tmp_formula_str = paste0('SOFR_IOER ~ Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + NetPosition_try + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins')
    tmp_formula_str = paste0('SOFR_IOER ~ Repo_balances1 + QuarterEnd + tdiff50_recv1 + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + Corporation_Income_Taxes_trillion + DlrDepUnins')
    commandline = paste0('use_labels(data2work, rq(formula(', tmp_formula_str, '), data = data2work,  tau = ', tau_tem, ') )' )
    tmp_model = eval(parse(text = commandline))
    coeff_vector_dealer2[i] = tmp_model$coefficients[2]
    coeff_vector_time2[i] = tmp_model$coefficients[4]
    reg_list2[[i]] = tmp_model
    Rsquare_list2[i] = 1 - tmp_model$rho/tmp_model0$rho
    
}

i=0
for (tau_tem in tau_vector){
    print(tau_tem)
    i = i + 1
    
    tmp_model = reg_list[[i]]
    
    robust_se[[i]] = summary.rq(tmp_model, se='ker')$coefficients[,'Std. Error']
    se_vector_dealer[i] = robust_se[[i]]["`dealer opening balances`"]
    se_vector_time[i] = robust_se[[i]]["`median time of receives`"]
    
    tmp_model = reg_list2[[i]]
    
    robust_se2[[i]] = summary.rq(tmp_model, se='ker')$coefficients[,'Std. Error']
    se_vector_dealer2[i] = robust_se2[[i]]["`dealer opening balances`"]
    se_vector_time2[i] = robust_se2[[i]]["`median time of receives`"]
}


se_vector_dealerBT = se_vector_dealer
se_vector_dealerBT2 = se_vector_dealer2
se_vector_timeBT = se_vector_time
se_vector_timeBT2 = se_vector_time2

i=0
for (tau_tem in tau_vector){
    print(tau_tem)
    i = i + 1
    
    tmp_model = reg_list[[i]]
    
    robust_se[[i]] = summary.rq(tmp_model, se='boot')$coefficients[,'Std. Error']
    se_vector_dealerBT[i] = robust_se[[i]]["`dealer opening balances`"]
    se_vector_timeBT[i] = robust_se[[i]]["`median time of receives`"]
    
    tmp_model = reg_list2[[i]]
    
    robust_se2[[i]] = summary.rq(tmp_model, se='boot')$coefficients[,'Std. Error']
    se_vector_dealerBT2[i] = robust_se2[[i]]["`dealer opening balances`"]
    se_vector_timeBT2[i] = robust_se2[[i]]["`median time of receives`"]
}


graphics.off()
Data2plot = tibble(tau_vector, coeff_vector_dealer, coeff_vector_time, Rsquare_list, se_vector_dealer, se_vector_time, se_vector_dealerBT, se_vector_timeBT) 
Data2plot2 = tibble(tau_vector, coeff_vector_dealer2, coeff_vector_time2, Rsquare_list2, se_vector_dealer2, se_vector_time2, se_vector_dealerBT2, se_vector_timeBT2) 




graph1=
    ggplot(Data2plot, aes(x = tau_vector)) +
    # 
    geom_line(aes(y = coeff_vector_dealer,    color= "dealeropeningbalances"   ),  linewidth=0.8) +
    # geom_line(aes(y = coeff_vector_time,     color= "mediantimeofreceives"    ),  linewidth=0.8) +
    # 
    scale_x_continuous(
        name = "percentile",
        # breaks = date_breaks("1 hour"),
        # breaks=date_breaks('6 months'),
        # labels = date_format("%b %Y", tz = "EST"),
        labels = percent,
        limits = c(0.5,1.0),
        n.breaks = 10,
        expand = c(0,0)
    ) +
    # the following specification should not be changed
    scale_y_continuous(name = "coefficient",
                       expand = c(0,0),   # Align plots on y axis
                       limits = c(-45,0),
                       n.breaks = 10
    ) +
    scale_color_manual(name = "Legend",
                       values=c( "dealeropeningbalances" = "blue", "mediantimeofreceives" = "red"),
                       labels = c( "dealeropeningbalances" = "dealer opening balances", "mediantimeofreceives" = "median time of receives")
    ) +
    theme(
        text = element_text(size = 25),       
        axis.text = element_text(colour = "black", size = rel(0.8)),
        axis.text.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit = "in")   ),
        axis.text.y = element_text(margin = margin(t=0, r=0.1, b=0,l=0, unit = "in")   ),
        # 
        panel.background = element_rect(fill = "white"),
        axis.line = element_line(colour = "black", 
                                 linewidth = 0.8, linetype = "solid"),
        axis.ticks = element_line(colour = "black"),
        axis.ticks.length = unit(-0.05,"in"),
        axis.title.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit='in') ),
        # 
        plot.margin = unit(c(0.5,0.5,0.2,0.6), "in"),  # top right bottom left
        panel.grid.major = element_line(colour = "grey80"),
        # panel.grid.minor = element_line(colour = "grey80"),
        legend.position = c(0.13, 0.95), 
        legend.title = element_blank()
    ) 
# ggplotly()
graph1
ggsave('dealeropeningbalancesQuantileCoefficientNoIV.pdf',graph1,width=16,height=9, units = "in" )



graph1=
    ggplot(Data2plot, aes(x = tau_vector)) +
    # 
    # geom_line(aes(y = coeff_vector_dealer,    color= "dealeropeningbalances"   ),  linewidth=0.8) +
    geom_line(aes(y = coeff_vector_time,     color= "mediantimeofreceives"    ),  linewidth=0.8) +
    # 
    scale_x_continuous(
        name = "percentile",
        # breaks = date_breaks("1 hour"),
        # breaks=date_breaks('6 months'),
        # labels = date_format("%b %Y", tz = "EST"),
        labels = percent,
        limits = c(0.5,1.0),
        n.breaks = 10,
        expand = c(0,0)
    ) +
    # the following specification should not be changed
    scale_y_continuous(name = "coefficient",
                       expand = c(0,0),   # Align plots on y axis
                       limits = c(0.085,0.135),
                       n.breaks = 10
    ) +
    scale_color_manual(name = "Legend",
                       values=c( "dealeropeningbalances" = "blue", "mediantimeofreceives" = "red"),
                       labels = c( "dealeropeningbalances" = "dealer opening balances", "mediantimeofreceives" = "median time of receives")
    ) +
    theme(
        text = element_text(size = 25),       
        axis.text = element_text(colour = "black", size = rel(0.8)),
        axis.text.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit = "in")   ),
        axis.text.y = element_text(margin = margin(t=0, r=0.1, b=0,l=0, unit = "in")   ),
        # 
        panel.background = element_rect(fill = "white"),
        axis.line = element_line(colour = "black", 
                                 linewidth = 0.8, linetype = "solid"),
        axis.ticks = element_line(colour = "black"),
        axis.ticks.length = unit(-0.05,"in"),
        axis.title.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit='in') ),
        # 
        plot.margin = unit(c(0.5,0.5,0.2,0.6), "in"),  # top right bottom left
        panel.grid.major = element_line(colour = "grey80"),
        # panel.grid.minor = element_line(colour = "grey80"),
        legend.position = c(0.13, 0.95), 
        legend.title = element_blank()
    ) 
# ggplotly()
graph1
ggsave('mediantimeofreceivesNOIV.pdf',graph1,width=16,height=9, units = "in" )





graph1=
    ggplot(Data2plot, aes(x = tau_vector)) +
    # 
    # geom_line(aes(y = coeff_vector_dealer,    color= "dealeropeningbalances"   ),  linewidth=0.8) +
    # geom_line(aes(y = coeff_vector_time,     color= "mediantimeofreceives"    ),  linewidth=0.8) +
    geom_line(aes(y = Rsquare_list,     color= "Rsquare_list"    ),  linewidth=0.8) +
    # 
    scale_x_continuous(
        name = "percentile",
        # breaks = date_breaks("1 hour"),
        # breaks=date_breaks('6 months'),
        # labels = date_format("%b %Y", tz = "EST"),
        labels = percent,
        limits = c(0.5,1.0),
        n.breaks = 10,
        expand = c(0,0)
    ) +
    # the following specification should not be changed
    scale_y_continuous(name = "coefficient",
                       expand = c(0,0),   # Align plots on y axis
                       limits = c(0.35,0.45),
                       n.breaks = 10
    ) +
    scale_color_manual(name = "Legend",
                       values=c( "Rsquare_list" = "black"),
                       labels = c( "Rsquare_list" = TeX("Sodu $R^2$ of the quantile regression"))
    ) +
    theme(
        text = element_text(size = 25),       
        axis.text = element_text(colour = "black", size = rel(0.8)),
        axis.text.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit = "in")   ),
        axis.text.y = element_text(margin = margin(t=0, r=0.1, b=0,l=0, unit = "in")   ),
        # 
        panel.background = element_rect(fill = "white"),
        axis.line = element_line(colour = "black", 
                                 linewidth = 0.8, linetype = "solid"),
        axis.ticks = element_line(colour = "black"),
        axis.ticks.length = unit(-0.05,"in"),
        axis.title.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit='in') ),
        # 
        plot.margin = unit(c(0.5,0.5,0.2,0.6), "in"),  # top right bottom left
        panel.grid.major = element_line(colour = "grey80"),
        # panel.grid.minor = element_line(colour = "grey80"),
        legend.position = c(0.19, 0.95), 
        legend.title = element_blank()
    ) 
# ggplotly()
graph1
ggsave('Rsquare_listNOIV.pdf',graph1,width=16,height=9, units = "in" )



# --------------Using another regression model

library("gridExtra")


graph1=
    ggplot(Data2plot2, aes(x = tau_vector)) +
    # 
    geom_ribbon(aes(ymin=coeff_vector_dealer2 - se_vector_dealer2 ,ymax=coeff_vector_dealer2 + se_vector_dealer2 , fill="dealeropeningbalances"), alpha=0.5) +
    geom_line(aes(y = coeff_vector_dealer2,    color= "dealeropeningbalances"   ),  linewidth=1.2) +
    # geom_line(aes(y = coeff_vector_time,     color= "mediantimeofreceives"    ),  linewidth=0.8) +
    # 
    scale_x_continuous(
        name = "percentile",
        # breaks = date_breaks("1 hour"),
        # breaks=date_breaks('6 months'),
        # labels = date_format("%b %Y", tz = "EST"),
        labels = percent,
        limits = c(0.65,1.0),
        n.breaks = 10,
        expand = c(0,0)
    ) +
    # the following specification should not be changed
    scale_y_continuous(name = "coefficient of dealer opening balances",
                       expand = c(0,0),   # Align plots on y axis
                       # limits = c(-45,0),
                       n.breaks = 10
    ) +
    scale_color_manual(name = "Legend",
                       values=c( "dealeropeningbalances" = "blue", "mediantimeofreceives" = "red"),
                       labels = c( "dealeropeningbalances" = "dealer opening balances", "mediantimeofreceives" = "median time of receives")
    ) +
    scale_fill_manual(name = "Legend",
                       values=c( "dealeropeningbalances" = "pink", "mediantimeofreceives" = "red"),
                       labels = c( "dealeropeningbalances" = "dealer opening balances", "mediantimeofreceives" = "median time of receives")
    ) +
    theme(
        text = element_text(size = 25),       
        axis.text = element_text(colour = "black", size = rel(0.8)),
        axis.text.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit = "in")   ),
        axis.text.y = element_text(margin = margin(t=0, r=0.1, b=0,l=0, unit = "in")   ),
        # 
        panel.background = element_rect(fill = "white"),
        axis.line = element_line(colour = "black", 
                                 linewidth = 0.8, linetype = "solid"),
        axis.ticks = element_line(colour = "black"),
        axis.ticks.length = unit(-0.05,"in"),
        axis.title.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit='in') ),
        # 
        plot.margin = unit(c(0.5,0.5,0.2,0.6), "in"),  # top right bottom left
        panel.grid.major = element_line(colour = "grey80"),
        # panel.grid.minor = element_line(colour = "grey80"),
        legend.position = c(0.33, 0.15), 
        legend.title = element_blank()
    ) 
graph1
ggsave('dealeropeningbalancesQuantileCoefficientNoIV2.pdf',graph1,width=16,height=9, units = "in" )



graph2=
    ggplot(Data2plot2, aes(x = tau_vector)) +
    # 
    geom_ribbon(aes(ymin=coeff_vector_time2 - se_vector_time2 ,ymax=coeff_vector_time2 + se_vector_time2 , fill="mediantimeofreceives"), alpha=0.5) +
    geom_line(aes(y = coeff_vector_time2,     color= "mediantimeofreceives"    ),  linewidth=1.2) +
    # 
    scale_x_continuous(
        name = "percentile",
        labels = percent,
        limits = c(0.65,1.0),
        n.breaks = 10,
        expand = c(0,0)
    ) +
    # the following specification should not be changed
    scale_y_continuous(name = "coefficient of median time of receives",
                       expand = c(0,0),   # Align plots on y axis
                       # limits = c(0.06,0.14),
                       n.breaks = 10
    ) +
    scale_color_manual(name = "Legend",
                       values=c( "dealeropeningbalances" = "blue", "mediantimeofreceives" = "red"),
                       labels = c( "dealeropeningbalances" = "dealer opening balances", "mediantimeofreceives" = "median time of receives")
    ) +
    scale_fill_manual(name = "Legend",
                      values=c( "dealeropeningbalances" = "pink", "mediantimeofreceives" = "pink"),
                      labels = c( "dealeropeningbalances" = "dealer opening balances", "mediantimeofreceives" = "median time of receives")
    ) +
    theme(
        text = element_text(size = 25),       
        axis.text = element_text(colour = "black", size = rel(0.8)),
        axis.text.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit = "in")   ),
        axis.text.y = element_text(margin = margin(t=0, r=0.1, b=0,l=0, unit = "in")   ),
        # 
        panel.background = element_rect(fill = "white"),
        axis.line = element_line(colour = "black", 
                                 linewidth = 0.8, linetype = "solid"),
        axis.ticks = element_line(colour = "black"),
        axis.ticks.length = unit(-0.05,"in"),
        axis.title.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit='in') ),
        # 
        plot.margin = unit(c(0.5,0.5,0.2,0.6), "in"),  # top right bottom left
        panel.grid.major = element_line(colour = "grey80"),
        # panel.grid.minor = element_line(colour = "grey80"),
        legend.position = c(0.33, 0.95), 
        legend.title = element_blank()
    ) 
# ggplotly()
graph2
ggsave('mediantimeofreceivesNOIV2.pdf',graph1,width=16,height=9, units = "in" )



# Figure 4
graph3 = grid.arrange(graph1, graph2, ncol= 2)
ggsave('QuantileNOIV2.pdf',graph3,width=16,height=9, units = "in" )



graph1=
    ggplot(Data2plot2, aes(x = tau_vector)) +
    # 
    geom_line(aes(y = Rsquare_list2,     color= "Rsquare_list"    ),  linewidth=0.8) +
    # 
    scale_x_continuous(
        name = "percentile",
        labels = percent,
        limits = c(0.5,1.0),
        n.breaks = 10,
        expand = c(0,0)
    ) +
    # the following specification should not be changed
    scale_y_continuous(name = "coefficient",
                       expand = c(0,0),   # Align plots on y axis
                       limits = c(0.38,0.475),
                       n.breaks = 10
    ) +
    scale_color_manual(name = "Legend",
                       values=c( "Rsquare_list" = "black"),
                       labels = c( "Rsquare_list" = TeX("Sodu $R^2$ of the quantile regression"))
    ) +
    theme(
        text = element_text(size = 25),       
        axis.text = element_text(colour = "black", size = rel(0.8)),
        axis.text.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit = "in")   ),
        axis.text.y = element_text(margin = margin(t=0, r=0.1, b=0,l=0, unit = "in")   ),
        # 
        panel.background = element_rect(fill = "white"),
        axis.line = element_line(colour = "black", 
                                 linewidth = 0.8, linetype = "solid"),
        axis.ticks = element_line(colour = "black"),
        axis.ticks.length = unit(-0.05,"in"),
        axis.title.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit='in') ),
        # 
        plot.margin = unit(c(0.5,0.5,0.2,0.6), "in"),  # top right bottom left
        panel.grid.major = element_line(colour = "grey80"),
        # panel.grid.minor = element_line(colour = "grey80"),
        legend.position = c(0.19, 0.95), 
        legend.title = element_blank()
    ) 
# ggplotly()
graph1
ggsave('Rsquare_listNOIV2.pdf',graph1,width=16,height=9, units = "in" )




# Median receive time and other bank reserves--------------------------------------------
graphics.off()
rm( list = setdiff( ls(), c(keepVariables, "") ))
cat("\014")  

data2work = data2work1 %>% mutate(critical_flag = if_else( (SOFR_IOER >25 | GCF_IOER >25 ), 1, 0 ))


tem =  data2work %>% filter(critical_flag == 1) %>% 
    top_n(20, wt = SOFR_IOER) %>% 
    mutate(largest_20_SOFR_IOER_spike = 1) %>% 
    select(Date, SOFR_IOER, largest_20_SOFR_IOER_spike)
tem1 = tem %>% top_n(10, wt = SOFR_IOER) %>% 
    mutate(largest_10_SOFR_IOER_spike = 1) %>% 
    select(Date, largest_10_SOFR_IOER_spike)
tem = left_join(tem, tem1, by='Date') %>% select(Date, largest_10_SOFR_IOER_spike, largest_20_SOFR_IOER_spike)

Data2plot = left_join(data2work, tem, by='Date') 


graph=
    ggplot(Data2plot, aes(x = AO_balances1*1000, y = tdiff50_recv1 ) ) +
    geom_point(size=1.5, color = 'blue') +
    geom_point(data = subset(Data2plot, largest_10_SOFR_IOER_spike == 1 ),  size=2, color = 'red') +
    geom_smooth(method='lm', formula= y~x, se = FALSE, color = 'black', size = 2) +
    geom_point(data = subset(Data2plot, Date == ymd( '2019-09-17', tz = 'EST' ) ), aes(x = AO_balances1*1000, y = tdiff50_recv1 ),  size=2.5, color = 'red') +
    labs(
        # x = TeX('Average daily dealer-bank opening balances ($ billions)'),
        x = ('Other large bank balances ($ billions)'),
        y = 'Time from sample mean (minutes)') +
    coord_cartesian( expand = TRUE ) +
    theme(text = element_text(size = 27),       
          axis.text = element_text(colour = "black", size = rel(0.8)),
          axis.text.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit = "in")   ),
          axis.text.y = element_text(margin = margin(t=0, r=0.1, b=0,l=0, unit = "in")   ),
          # 
          panel.background = element_rect(fill = "white"),
          axis.line = element_line(colour = "black", 
                                   size = 0.1, linetype = "solid"),
          axis.ticks = element_line(colour = "black"),
          axis.ticks.length = unit(-0.05,"in"),
          axis.title.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit='in') ),
          # 
          plot.margin = unit(c(0.5,0.5,0.2,0.6), "in"),  # top right bottom left
          panel.grid.major = element_line(colour = "grey80"),
          # panel.grid.minor = element_line(colour = "grey80")    
          # legend.title =element_blank(), 
          # legend.position= c(0.2, 0.95)
    ) 

graph
ggsave('median_payment_time_vs_other_banks_balances_top10_SOFR_SpikeDays.pdf',graph,width=16,height=9, units = "in" )

# Figure 3
graph=
    ggplot(Data2plot, aes(x = AO_balances1*1000, y = tdiff50_recv1 ) ) +
    geom_point(size=1.5, color = 'blue') +
    geom_point(data = subset(Data2plot, largest_20_SOFR_IOER_spike == 1 ),  size=2, color = 'red') +
    geom_smooth(method='lm', formula= y~x, se = FALSE, color = 'black', size = 2) +
    geom_point(data = subset(Data2plot, Date == ymd( '2019-09-17', tz = 'EST' ) ), aes(x = AO_balances1*1000, y = tdiff50_recv1 ),  size=2.5, color = 'red') +
    labs(
        # x = TeX('Average daily dealer-bank opening balances ($ billions)'),
        x = ('Other large bank balances ($ billions)'),
        # y = 'Time from sample mean (minutes)'
        y = 'Payment delay (minutes)'
        ) +
    # scale_x_continuous(  expand = c(0,0)    ) +
    # scale_y_continuous(
    #     # the following specification should not be changed
    #     expand = c(0,0),   # Align plots on y axis
    # ) +    
    coord_cartesian( expand = TRUE ) +
    # coord_cartesian(xlim = c(min(x)-10,max(x)+10), ylim = c(min(y)-10,max(y)+10) ) +
    # labs(labels = c("Marketable treasuries outstanding", "projected")) +
    theme(text = element_text(size = 27),       
          axis.text = element_text(colour = "black", size = rel(0.8)),
          axis.text.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit = "in")   ),
          axis.text.y = element_text(margin = margin(t=0, r=0.1, b=0,l=0, unit = "in")   ),
          # 
          panel.background = element_rect(fill = "white"),
          axis.line = element_line(colour = "black", 
                                   size = 0.1, linetype = "solid"),
          axis.ticks = element_line(colour = "black"),
          axis.ticks.length = unit(-0.05,"in"),
          axis.title.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit='in') ),
          # 
          plot.margin = unit(c(0.5,0.5,0.2,0.6), "in"),  # top right bottom left
          panel.grid.major = element_line(colour = "grey80"),
          # panel.grid.minor = element_line(colour = "grey80")    
          # legend.title =element_blank(), 
          # legend.position= c(0.2, 0.95)
    ) 

graph
ggsave('median_payment_time_vs_other_banks_balances_top20_SOFR_SpikeDays2024.pdf',graph,width=16,height=9, units = "in" )




# ---------------------Time series of reserves------------------------------------------------------
# Figure 2
rm( list = setdiff( ls(), c(keepVariables, "") ))
Data2plot = data2work0
tem1 = read.csv(file = 'TotalReserves2015_2020Sept.csv')  # in millions of dollars
tem1 %<>% mutate(TotalReserves = RESBALNSW/1000000,  # to trillions of dollars   
                Date = mdy(Date, tz = 'EST') ) %>% 
    select(Date, TotalReserves)



tem2 = read.csv(file = 'TotalReserveMonthly2015to2023.csv')  # in billions of dollars
tem2 %<>% mutate(TotalReserves = TOTRESNS/1000,  # to trillions of dollars   
                Date = ymd(DATE, tz = 'EST') ) %>% 
    select(Date, TotalReserves) %>%
    filter(Date >= ymd('2015-01-01') )

tem2[1,]$Date = ymd('2015-01-02', tz = 'EST')
tem3 = rbind(tem2[1,], tem1 )

tem2 %<>% filter(Date > tem3[nrow(tem3),1] )

tem = rbind(tem3, tem2)
Data2plot2tem0 = data2work0 %>% 
    arrange(Date) %>% 
    mutate(DealerReserves =   Repo_balances1,
           DealerOtherReserves = AO_balances1 + Repo_balances1, 
           # TotalReserves = TotalBalBil, 
    ) %>%
    select(Date, DealerReserves, DealerOtherReserves)


Data2plot2tem = full_join(Data2plot2tem0, tem, by = 'Date') %>% arrange(Date) %>% filter(Date <= ymd('2023-03-13') )
# Data2plot2tem$TotalReserves[1] = 2.587691


Data2plot = Data2plot2tem %>% 
    arrange(Date) %>%
    fill(DealerReserves, DealerOtherReserves, TotalReserves, .direction = "down") %>% 
    mutate(DealerReserves = DealerReserves * 1000, 
           DealerOtherReserves = DealerOtherReserves * 1000,
           TotalReserves = TotalReserves * 1000)


graph1=
    ggplot(Data2plot, aes(x = Date)) +
    geom_line(aes(y = DealerReserves), color='blue',size=0.3) +
    geom_line(aes(y = DealerOtherReserves), color='red',size=0.3) +
    geom_line(aes(y = TotalReserves), color='green',size=0.3) +
    geom_ribbon(aes(ymin=0,ymax=DealerReserves, fill="large dealer-banks"), alpha=0.8) +
    geom_ribbon(aes(ymin=DealerReserves,ymax=DealerOtherReserves, fill="other large banks"), alpha=0.6) +
    geom_ribbon(aes(ymin=DealerOtherReserves,ymax=TotalReserves, fill="other"), alpha=0.7) +
    scale_x_datetime(name = element_blank(),
                     # breaks = date_breaks("1 hour"),
                     breaks=date_breaks('6 months'),
                     minor_breaks = date_breaks('6 months'),
                     labels = date_format("%m/%Y", tz = "EST"),
                     # limits = c(Regression4$Date[2],Regression4$ Date[328])
                     expand = c(0,0)
    ) +
    # the following specification should not be changed
    scale_y_continuous(name = "Reserve balances ($ billions)",
                       expand = c(0,0),   # Align plots on y axis
                       # limits = c(0,27),
                       # n.breaks = 6
    ) +
    scale_fill_manual(values=c( "large dealer-banks" = "blue", "other large banks" = "red", "other" = "green"),
                      breaks=c("large dealer-banks", "other large banks" , "other" )
    ) +
    theme(
        text = element_text(size = 25),
        axis.text = element_text(colour = "black", size = rel(0.8)),
        axis.text.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit = "in"), angle=45, hjust = 1  ),
        axis.text.y = element_text(margin = margin(t=0, r=0.1, b=0,l=0, unit = "in")   ),
        #
        panel.background = element_rect(fill = "white"),
        axis.line = element_line(colour = "black",
                                 size = 0.8, linetype = "solid"),
        axis.ticks = element_line(colour = "black"),
        axis.ticks.length = unit(-0.05,"in"),
        axis.title.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit='in') ),
        #
        plot.margin = unit(c(0.5,0.5,0.2,0.6), "in"),  # top right bottom left
        # panel.grid.major = element_line(colour = "grey80"),
        # panel.grid.minor = element_line(colour = "grey80")
        legend.position = c(0.15, 0.95),
        legend.title = element_blank()
    )
# ggplotly()
graph1
ggsave('ReserveBreakdown2015toMarch2023.pdf',graph1,width=16,height=9, units = "in" )


# 





# ----------Figure 11----------------------------------------
# using Kaplan Sun mathed

tau_vector = c(
    51	,
    53	,
    55	,
    57	,
    59	,
    61	,
    63	,
    65	,
    67	,
    69	,
    71	,
    73	,
    75	,
    77	,
    79	,
    81	,
    83	,
    85	,
    # 89	,
    91	,
    93	,
    95	,
    97	,
    99	
)

# The following data are from the Stata File in the folder Table2Table14

alpha_vec = c(
    -27.38639	,
    -27.46456	,
    -27.54624	,
    -29.17885	,
    -27.03563	,
    -26.90972	,
    -26.54742	,
    -25.92859	,
    -25.21829	,
    -24.67587	,
    -24.06095	,
    -22.87164	,
    -22.71825	,
    -22.83229	,
    -22.24422	,
    -21.79577	,
    -21.56462	,
    -21.29877	,
    -26.35772	,
    -31.04483	,
    -35.37141	,
    -44.15714	,
    -74.11317	
)

std_vector_dealer = c(
    4.172786	,
    3.997233	,
    3.879954	,
    4.067715	,
    3.811937	,
    3.887682	,
    3.935058	,
    3.979571	,
    4.114883	,
    4.307986	,
    4.603999	,
    5.226453	,
    5.776907	,
    6.511899	,
    6.773187	,
    7.418922	,
    8.112322	,
    9.460415	,
    11.88593	,
    10.68988	,
    11.56141	,
    10.12188	,
    31.34177	
    )
 
Data2plot = tibble(tau_vector=tau_vector/100, coeff_vector_dealer = alpha_vec, std_vector_dealer, Rsquare_list = 0) 
# Data2plot %<>% filter(tau_vector>=0.6)
graph1=
    ggplot(Data2plot, aes(x = tau_vector)) +
    # 
    geom_line(aes(y = coeff_vector_dealer,    color= "dealeropeningbalances"   ),  linewidth=0.8) +
    geom_ribbon(aes(ymin=coeff_vector_dealer-std_vector_dealer,ymax=coeff_vector_dealer+std_vector_dealer, fill="large dealer-banks"), alpha=0.3) +
    # 
    scale_x_continuous(
        name = "Sensitivity to predicted dealer opening balances",
        # breaks = date_breaks("1 hour"),
        # breaks=date_breaks('6 months'),
        # labels = date_format("%b %Y", tz = "EST"),
        labels = percent,
        # limits = c(50,100),
        n.breaks = 10,
        expand = c(0,0)
    ) +
    # the following specification should not be changed
    scale_y_continuous(name = "coefficient",
                       expand = c(0,0),   # Align plots on y axis
                       # limits = c(-80,-15),
                       n.breaks = 10
    ) +
    scale_color_manual(name = "Legend",
                       values=c( "dealeropeningbalances" = "blue", "mediantimeofreceives" = "red"),
                       labels = c( "dealeropeningbalances" = "predicted dealer opening balances", "mediantimeofreceives" = "median time of receives")
    ) +
    # scale_color_discrete(name = "Legend") +
    # scale_linetype_manual(name = "Legend",
    #                       values=c( "SOFR_IOER" = "solid", "GCF_IOER" = "longdash", "TGCR_IOER" = "dashed"),
    #                       labels = c( "SOFR_IOER" = "SOFR-IOER", "GCF_IOER" = "GCF-IOER", "TGCR_IOER" = "TGCR-IOER" )
    #                       ) +
    theme(
        text = element_text(size = 25),       
        axis.text = element_text(colour = "black", size = rel(0.8)),
        axis.text.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit = "in")   ),
        axis.text.y = element_text(margin = margin(t=0, r=0.1, b=0,l=0, unit = "in")   ),
        # 
        panel.background = element_rect(fill = "white"),
        axis.line = element_line(colour = "black", 
                                 linewidth = 0.8, linetype = "solid"),
        axis.ticks = element_line(colour = "black"),
        axis.ticks.length = unit(-0.05,"in"),
        axis.title.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit='in') ),
        # 
        plot.margin = unit(c(0.5,0.5,0.2,0.6), "in"),  # top right bottom left
        panel.grid.major = element_line(colour = "grey80"),
        # panel.grid.minor = element_line(colour = "grey80"),
        # legend.position = c(0.2, 0.15), 
        # legend.position = element_blank(), 
        legend.position="none",
        legend.title = element_blank()
    ) 
# ggplotly()
graph1
ggsave('dealeropeningbalancesQuantileCoefficientIV202406.pdf',graph1,width=16,height=9, units = "in" )










# ---------------------Main Regression with lagged payment timing -----------------------------------------------------------------------
# First definition
# Table 3
rm( list = setdiff( ls(), c(keepVariables, "") ))
data2work = data2work1 %<>% 
    mutate(Ave_tdiff50_recv1_3 = zoo::rollmean(tdiff50_recv1, k = 3, align = "right", fill = NA),
           Ave_tdiff50_recv1_5 = zoo::rollmean(tdiff50_recv1, k = 5, align = "right", fill = NA),
           Ave_tdiff50_recv1_7 = zoo::rollmean(tdiff50_recv1, k = 7, align = "right", fill = NA),
           Ave_tdiff50_recv1_9 = zoo::rollmean(tdiff50_recv1, k = 9, align = "right", fill = NA),
           Ave_tdiff50_recv1_11 = zoo::rollmean(tdiff50_recv1, k = 11, align = "right", fill = NA),
           Ave_tdiff50_recv1_13 = zoo::rollmean(tdiff50_recv1, k = 13, align = "right", fill = NA),
           Ave_tdiff50_recv1_15 = zoo::rollmean(tdiff50_recv1, k = 15, align = "right", fill = NA),
           Ave_tdiff50_recv1_17 = zoo::rollmean(tdiff50_recv1, k = 17, align = "right", fill = NA),
           Ave_tdiff50_recv1_19 = zoo::rollmean(tdiff50_recv1, k = 19, align = "right", fill = NA),
           Ave_tdiff50_recv1_21 = zoo::rollmean(tdiff50_recv1, k = 21, align = "right", fill = NA),
           AveBalancesPreviousWeeklag3 = lag(Ave_tdiff50_recv1_3, 2),
           AveBalancesPreviousWeeklag5 = lag(Ave_tdiff50_recv1_5, 2),
           AveBalancesPreviousWeeklag7 = lag(Ave_tdiff50_recv1_7, 2),
           AveBalancesPreviousWeeklag9 = lag(Ave_tdiff50_recv1_9, 2),
           AveBalancesPreviousWeeklag11 = lag(Ave_tdiff50_recv1_11, 2),
           AveBalancesPreviousWeeklag13 = lag(Ave_tdiff50_recv1_13, 2),
           AveBalancesPreviousWeeklag15 = lag(Ave_tdiff50_recv1_15, 2),
           AveBalancesPreviousWeeklag17 = lag(Ave_tdiff50_recv1_17, 2),
           AveBalancesPreviousWeeklag19 = lag(Ave_tdiff50_recv1_19, 2),
           AveBalancesPreviousWeeklag21 = lag(Ave_tdiff50_recv1_21, 2),
           
           lag2_Repo_balances1 = lag(Repo_balances1, 2)
           )



data2work =
    apply_labels(
        data2work,
        AveBalancesPreviousWeeklag9 = "Average Payment Delay t-10 to t-2",
        lag2_Repo_balances1 = "dealer opening balances t-2"
    )


regressor_list = c(
    # 'lag2_Repo_balances1 + QuarterEnd',
    'lag2_Repo_balances1 + AveBalancesPreviousWeeklag9 + QuarterEnd',
    'AveBalancesPreviousWeeklag9 + QuarterEnd ',
    'lag2_Repo_balances1 + QuarterEnd + AveBalancesPreviousWeeklag9 + Tbills_outstanding + DlrDepUnins',
    'lag2_Repo_balances1 + QuarterEnd + AveBalancesPreviousWeeklag9 + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins + Corporation_Income_Taxes_trillion'
)


dep_list = c('SOFR_IOER')

reg_list = vector(mode = 'list', length = length(regressor_list) * length(dep_list) )
robust_se = vector(mode = 'list', length = length(regressor_list) * length(dep_list)  )
x_vars = c()
i = 1

for (a1 in regressor_list) {
    for (a2 in dep_list) {
        tmp_formula_str = paste0(a2, ' ~ ', a1)
        commandline = paste0('use_labels(data2work, lm(formula(', tmp_formula_str, '), data = data2work ) )' )
        tmp_model = eval(parse(text = commandline))
        # tmp_model = lm(as.formula(tmp_formula_str), data = data2work )
        # tmp_model = lm(as.formula(tmp_formula_str), data = Data2regress)
        cov1         = vcovHC(tmp_model, type = "HC1")
        robust_se[[i]] = sqrt(diag(cov1))
        reg_list[[i]] = tmp_model
        x_vars = c(x_vars, a2)
        i = i + 1
    }
}
note.latex = paste0(
    "\\multicolumn{", 
    length(reg_list)+1,
    "}{l}{{Note:} Standard errors are adjusted for heteroskedasticity. $^{*}p<0.1$; $^{**}p<0.05$; $^{***}p<0.01$.} \\\\ \\multicolumn{",
    length(reg_list)+1,
    "}{l}{ \\qquad \\quad Constant included for each specification.}")

res = remove_backticks( 
    stargazer(
        reg_list, 
        title="Results", 
        align=TRUE, 
        omit.stat = c("f"), 
        dep.var.caption  =  paste0("Dependent variable: ", eval(parse(text = paste0('var_lab(data2work$',dep_list[1],')'))) ),
        dep.var.labels.include = FALSE,
        df = FALSE,
        digits = 5,
        # notes        = TeX("Sometimes you just have to start over."), 
        # notes.append = FALSE,
        se = robust_se))
res[grepl("Note",res)] = note.latex
res[grepl("\\label\\{\\}",res)] = "%"
res[grepl("\\caption\\{",res)] = "%"

res[grepl("Table created by stargazer",res)] = ""
res[grepl("\\% Date and time:",res)] = ""
res[grepl("\\% Requires LaTeX packages: dcolumn",res)] = ""

regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res)) =  lapply(regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res) ) , mysignif )
# 
res[grepl("begin\\{table\\}",res)] = begintable.latex
res[grepl("end\\{table\\}",res)] = endtable.latex
res = gsub("R\\$\\^\\{2\\}\\$", "\\$R\\^\\{2\\}\\$", res)
cat("\014")  
write_clip(res)
writeLines(res)


# ----------------------Quantile Regression ---Average Payment Delay-------------------------
# Table 9
library(quantreg)
rm( list = setdiff( ls(), c(keepVariables, "") ))
data2work = data2work1 %<>% 
    mutate(Ave_tdiff50_recv1_3 = zoo::rollmean(tdiff50_recv1, k = 3, align = "right", fill = NA),
           Ave_tdiff50_recv1_5 = zoo::rollmean(tdiff50_recv1, k = 5, align = "right", fill = NA),
           Ave_tdiff50_recv1_7 = zoo::rollmean(tdiff50_recv1, k = 7, align = "right", fill = NA),
           Ave_tdiff50_recv1_9 = zoo::rollmean(tdiff50_recv1, k = 9, align = "right", fill = NA),
           Ave_tdiff50_recv1_11 = zoo::rollmean(tdiff50_recv1, k = 11, align = "right", fill = NA),
           Ave_tdiff50_recv1_13 = zoo::rollmean(tdiff50_recv1, k = 13, align = "right", fill = NA),
           Ave_tdiff50_recv1_15 = zoo::rollmean(tdiff50_recv1, k = 15, align = "right", fill = NA),
           Ave_tdiff50_recv1_17 = zoo::rollmean(tdiff50_recv1, k = 17, align = "right", fill = NA),
           Ave_tdiff50_recv1_19 = zoo::rollmean(tdiff50_recv1, k = 19, align = "right", fill = NA),
           Ave_tdiff50_recv1_21 = zoo::rollmean(tdiff50_recv1, k = 21, align = "right", fill = NA),
           AveBalancesPreviousWeeklag3 = lag(Ave_tdiff50_recv1_3, 2),
           AveBalancesPreviousWeeklag5 = lag(Ave_tdiff50_recv1_5, 2),
           AveBalancesPreviousWeeklag7 = lag(Ave_tdiff50_recv1_7, 2),
           AveBalancesPreviousWeeklag9 = lag(Ave_tdiff50_recv1_9, 2),
           AveBalancesPreviousWeeklag11 = lag(Ave_tdiff50_recv1_11, 2),
           AveBalancesPreviousWeeklag13 = lag(Ave_tdiff50_recv1_13, 2),
           AveBalancesPreviousWeeklag15 = lag(Ave_tdiff50_recv1_15, 2),
           AveBalancesPreviousWeeklag17 = lag(Ave_tdiff50_recv1_17, 2),
           AveBalancesPreviousWeeklag19 = lag(Ave_tdiff50_recv1_19, 2),
           AveBalancesPreviousWeeklag21 = lag(Ave_tdiff50_recv1_21, 2),
           
           lag2_Repo_balances1 = lag(Repo_balances1, 2)
    )

data2work =
    apply_labels(
        data2work,
        AveBalancesPreviousWeeklag9 = "Average Payment Delay t-10 to t-2",
        lag2_Repo_balances1 = "dealer opening balances t-2"
    )

tem = data2work %>% summarize(
    AveBalancesPreviousWeeklag3_mean=mean(AveBalancesPreviousWeeklag3, na.rm = TRUE),
    AveBalancesPreviousWeeklag5_mean=mean(AveBalancesPreviousWeeklag5, na.rm = TRUE),
    AveBalancesPreviousWeeklag7_mean=mean(AveBalancesPreviousWeeklag7, na.rm = TRUE),
    AveBalancesPreviousWeeklag9_mean=mean(AveBalancesPreviousWeeklag9, na.rm = TRUE),
    AveBalancesPreviousWeeklag11_mean=mean(AveBalancesPreviousWeeklag11, na.rm = TRUE),
    AveBalancesPreviousWeeklag13_mean=mean(AveBalancesPreviousWeeklag13, na.rm = TRUE),
    AveBalancesPreviousWeeklag15_mean=mean(AveBalancesPreviousWeeklag15, na.rm = TRUE),
    AveBalancesPreviousWeeklag17_mean=mean(AveBalancesPreviousWeeklag17, na.rm = TRUE),
    AveBalancesPreviousWeeklag19_mean=mean(AveBalancesPreviousWeeklag19, na.rm = TRUE),
    AveBalancesPreviousWeeklag21_mean=mean(AveBalancesPreviousWeeklag21, na.rm = TRUE),
    
    AveBalancesPreviousWeeklag3_sd=sd(AveBalancesPreviousWeeklag3, na.rm = TRUE),
    AveBalancesPreviousWeeklag5_sd=sd(AveBalancesPreviousWeeklag5, na.rm = TRUE),
    AveBalancesPreviousWeeklag7_sd=sd(AveBalancesPreviousWeeklag7, na.rm = TRUE),
    AveBalancesPreviousWeeklag9_sd=sd(AveBalancesPreviousWeeklag9, na.rm = TRUE),
    AveBalancesPreviousWeeklag11_sd=sd(AveBalancesPreviousWeeklag11, na.rm = TRUE),
    AveBalancesPreviousWeeklag13_sd=sd(AveBalancesPreviousWeeklag13, na.rm = TRUE),
    AveBalancesPreviousWeeklag15_sd=sd(AveBalancesPreviousWeeklag15, na.rm = TRUE),
    AveBalancesPreviousWeeklag17_sd=sd(AveBalancesPreviousWeeklag17, na.rm = TRUE),
    AveBalancesPreviousWeeklag19_sd=sd(AveBalancesPreviousWeeklag19, na.rm = TRUE),
    AveBalancesPreviousWeeklag21_sd=sd(AveBalancesPreviousWeeklag21, na.rm = TRUE),
    
    AveBalancesPreviousWeeklag3_min=min(AveBalancesPreviousWeeklag3, na.rm = TRUE),
    AveBalancesPreviousWeeklag5_min=min(AveBalancesPreviousWeeklag5, na.rm = TRUE),
    AveBalancesPreviousWeeklag7_min=min(AveBalancesPreviousWeeklag7, na.rm = TRUE),
    AveBalancesPreviousWeeklag9_min=min(AveBalancesPreviousWeeklag9, na.rm = TRUE),
    AveBalancesPreviousWeeklag11_min=min(AveBalancesPreviousWeeklag11, na.rm = TRUE),
    AveBalancesPreviousWeeklag13_min=min(AveBalancesPreviousWeeklag13, na.rm = TRUE),
    AveBalancesPreviousWeeklag15_min=min(AveBalancesPreviousWeeklag15, na.rm = TRUE),
    AveBalancesPreviousWeeklag17_min=min(AveBalancesPreviousWeeklag17, na.rm = TRUE),
    AveBalancesPreviousWeeklag19_min=min(AveBalancesPreviousWeeklag19, na.rm = TRUE),
    AveBalancesPreviousWeeklag21_min=min(AveBalancesPreviousWeeklag21, na.rm = TRUE),
    
    AveBalancesPreviousWeeklag3_max=max(AveBalancesPreviousWeeklag3, na.rm = TRUE),
    AveBalancesPreviousWeeklag5_max=max(AveBalancesPreviousWeeklag5, na.rm = TRUE),
    AveBalancesPreviousWeeklag7_max=max(AveBalancesPreviousWeeklag7, na.rm = TRUE),
    AveBalancesPreviousWeeklag9_max=max(AveBalancesPreviousWeeklag9, na.rm = TRUE),
    AveBalancesPreviousWeeklag11_max=max(AveBalancesPreviousWeeklag11, na.rm = TRUE),
    AveBalancesPreviousWeeklag13_max=max(AveBalancesPreviousWeeklag13, na.rm = TRUE),
    AveBalancesPreviousWeeklag15_max=max(AveBalancesPreviousWeeklag15, na.rm = TRUE),
    AveBalancesPreviousWeeklag17_max=max(AveBalancesPreviousWeeklag17, na.rm = TRUE),
    AveBalancesPreviousWeeklag19_max=max(AveBalancesPreviousWeeklag19, na.rm = TRUE),
    AveBalancesPreviousWeeklag21_max=max(AveBalancesPreviousWeeklag21, na.rm = TRUE)
)


regressor_list = c(
    'lag2_Repo_balances1 + AveBalancesPreviousWeeklag9 + QuarterEnd',
    'QuarterEnd + AveBalancesPreviousWeeklag9',
    'lag2_Repo_balances1 + QuarterEnd + AveBalancesPreviousWeeklag9 + Tbills_outstanding + DlrDepUnins',
    'lag2_Repo_balances1 + QuarterEnd + AveBalancesPreviousWeeklag9 + Tbills_outstanding + total_redemptions_all + total_issuance_bills + coupon_issuances + DlrDepUnins + Corporation_Income_Taxes_trillion'
)



mean(data2work$AveBalancesPreviousWeeklag9, na.rm = TRUE)
sd(data2work$AveBalancesPreviousWeeklag9, na.rm = TRUE)
sd(data2work$tdiff50_recv1, na.rm = TRUE)
min(data2work$AveBalancesPreviousWeeklag9, na.rm = TRUE)
max(data2work$AveBalancesPreviousWeeklag9, na.rm = TRUE)


dep_list = c('SOFR_IOER')

reg_list = vector(mode = 'list', length = length(regressor_list) * length(dep_list) )
robust_se = vector(mode = 'list', length = length(regressor_list) * length(dep_list) )
Rsquare_list = vector(mode = 'character', length = length(regressor_list) * length(dep_list) )

tmp_model0 = rq(SOFR_IOER ~ 1, data = data2work,  tau = 0.99)


x_vars = c()
i = 1
for (a1 in regressor_list) {
    for (a2 in dep_list) {
        tmp_formula_str = paste0(a2, ' ~ ', a1)
        commandline = paste0('use_labels(data2work, rq(formula(', tmp_formula_str, '), data = data2work,  tau = 0.99) )' )
        tmp_model = eval(parse(text = commandline))
        reg_list[[i]] = tmp_model
        # cov1         = vcovHC(tmp_model, type = "HC1")
        # robust_se[[i]] = sqrt(diag(cov1))
        Rsquare_list[[i]] = 1 - tmp_model$rho/tmp_model0$rho
        x_vars = c(x_vars, a2)
        i = i + 1
    }
}

note.latex = paste0(
    "\\multicolumn{",
    length(reg_list)+1,
    "}{l}{{Note:} Constant included for each specification.}")


res = remove_backticks(
    stargazer(
        reg_list,
        title="Results",
        align=TRUE,
        omit.stat = c("f"),
        dep.var.caption  =  paste0("Dependent variable: ", eval(parse(text = paste0('var_lab(data2work$',dep_list[1],')'))) ),
        dep.var.labels.include = FALSE,
        df = FALSE,
        digits = 5,
        # notes        = TeX("Sometimes you just have to start over."),
        # notes.append = FALSE,
        # se = robust_se
        add.lines=list( c("pseudo-$R^{1}$", Rsquare_list) ), 
        table.layout = "=lc#-t-sa-n"
    )
)
res[grepl("Note",res)] = note.latex
res[grepl("\\label\\{\\}",res)] = "%"
res[grepl("\\caption\\{",res)] = "%"
res[grepl("Table created by stargazer v.5.2.3 by Marek Hlavac",res)] = ""
res[grepl("\\% Date and time:",res)] = ""
res[grepl("\\% Requires LaTeX packages: dcolumn",res)] = ""
regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res)) =  lapply(regmatches(res, gregexpr("(-?)(\\d+(\\,\\d+)?(\\.\\d+)?)", res) ) , mysignif )
res[grepl("begin\\{table\\}",res)] = begintable.latex
res[grepl("end\\{table\\}",res)] = endtable.latex
cat("\014")  
write_clip(res)
writeLines(res)
# End








# ---------------------Time series plot of spikes------------------------------------------------------
graphics.off()
Data2plot = data2work0 %>% filter( !is.na(SOFR_IOER) & !is.na(GCF_IOER) ) %>% 
    mutate(SOFR_IOER_trancated = pmin(SOFR_IOER, 200), 
           GCF_IOER_trancated  = pmin(GCF_IOER,  200),
           TGCR_IOER_trancated = pmin(TGCR_IOER, 200))

# Figure 10
graphics.off()
Data2plot_tem = Data2plot %>% filter(Date > ymd('2018-01-01', tz = "EST") )%>% filter(Date < ymd('2020-9-5', tz = "EST") )
graph1=
    ggplot(Data2plot_tem, aes(x = Date)) +
    # 
    geom_line(aes(y = SOFR_IOER_trancated,    color= "SOFR_IOER_trancated"   ),  size=0.8) +
    geom_line(aes(y = GCF_IOER_trancated,     color= "GCF_IOER_trancated"    ),  size=0.8) +
    # 
    scale_x_datetime(
        # name = element_blank(), 
        # breaks = date_breaks("1 hour"),
        breaks=date_breaks('6 months'),
        labels = date_format("%b %Y", tz = "EST"),
        # limits = c(Regression4$Date[2],Regression4$ Date[328])
        expand = c(0,0)
    ) +
    # the following specification should not be changed
    scale_y_continuous(name = "Spread (basis points)",
                       expand = c(0,0),   # Align plots on y axis
                       # limits = c(0,27),
                       # n.breaks = 6
    ) +
    scale_color_manual(name = "Legend",
                       values=c( "SOFR_IOER_trancated" = "blue", "GCF_IOER_trancated" = "red"),
                       labels = c( "SOFR_IOER_trancated" = "SOFR-IOR (truncated)", "GCF_IOER_trancated" = "GCF-IOR (truncated)")
    ) +
    # scale_color_discrete(name = "Legend") +
    # scale_linetype_manual(name = "Legend",
    #                       values=c( "SOFR_IOER" = "solid", "GCF_IOER" = "longdash", "TGCR_IOER" = "dashed"),
    #                       labels = c( "SOFR_IOER" = "SOFR-IOER", "GCF_IOER" = "GCF-IOER", "TGCR_IOER" = "TGCR-IOER" )
    #                       ) +
    theme(
        text = element_text(size = 25),       
        axis.text = element_text(colour = "black", size = rel(0.8)),
        axis.text.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit = "in")   ),
        axis.text.y = element_text(margin = margin(t=0, r=0.1, b=0,l=0, unit = "in")   ),
        # 
        panel.background = element_rect(fill = "white"),
        axis.line = element_line(colour = "black", 
                                 size = 0.8, linetype = "solid"),
        axis.ticks = element_line(colour = "black"),
        axis.ticks.length = unit(-0.05,"in"),
        axis.title.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit='in') ),
        # 
        plot.margin = unit(c(0.5,0.5,0.2,0.6), "in"),  # top right bottom left
        panel.grid.major = element_line(colour = "grey80"),
        # panel.grid.minor = element_line(colour = "grey80"),
        legend.position = c(0.13, 0.95), 
        legend.title = element_blank()
    ) 
# ggplotly()
graph1
# ggsave('SOFR_IOERvsGCF_IOER_2018_truncated.pdf',graph1,width=16,height=9, units = "in" )






# Appendix Figure 9-----------------------------------------
rm( list = setdiff( ls(), c(keepVariables, "") ))
cat("\014")  
graphics.off()
dayoverdraft_raw1 = read_excel("DaylightPeakOverdraft.xlsx", sheet = "Data2010to2023", skip = 0 )
dayoverdraft_raw2 = read_excel("DaylightPeakOverdraft.xlsx", sheet = "Data2000to2009", skip = 0 )
dayoverdraft1 = dayoverdraft_raw1 %>% 
    mutate(EndingDate = Reserve_Maintenance_Period_Ending,  
           OverdraftFund = `Intraday Peak Overdrafts2_(million $)__Total`/1000, # from million to billion 
           EndingDate = ymd(EndingDate, tz = 'EST') ) %>% 
    select(EndingDate, OverdraftFund)
dayoverdraft2 = dayoverdraft_raw2 %>% 
    mutate(EndingDate = Reserve_Maintenance_Period_Ending,  
           OverdraftFund = `Intraday Peak Overdrafts_(million $)_Total`/1000, # from million to billion 
           EndingDate = ymd(EndingDate, tz = 'EST') ) %>% 
    select(EndingDate, OverdraftFund)
dayoverdraft = rbind(dayoverdraft1, dayoverdraft2) 
dayoverdraft %<>% mutate( EndingDate = ymd(EndingDate, tz = 'EST'),
                          Date =  EndingDate ) 
# rm( list = setdiff( ls(), c("remove_backticks", "dayoverdraft" )  ))
cat("\014")  
rm(dayoverdraft_raw1, dayoverdraft_raw2, dayoverdraft1, dayoverdraft2)

tem1 = dayoverdraft %>% select(EndingDate, OverdraftFund) %>% filter(EndingDate <= ymd('2009-12-31', tz = 'EST') )

tem2 = dayoverdraft %>% select(EndingDate, OverdraftFund) %>% filter(EndingDate >= ymd('2009-12-31', tz = 'EST') )

mean(tem2$OverdraftFund)

Data2plot = data2work0
Data2plot_tem = Data2plot %>% 
    mutate( Total_balances1 = (AO_balances1 + Repo_balances1)* 1000    ) %>% 
    select(Date, Total_balances1)

dayoverdraft4plot1 = left_join(Data2plot_tem, dayoverdraft, by = 'Date') %>% 
    fill(EndingDate, .direction = "up") %>%
    # fill(EndingDate, .direction = "down") %>% 
    group_by(EndingDate) %>% 
    mutate(averagedealerbalance = mean(Total_balances1, trim = 0, na.rm = TRUE)  ) %>% 
    select(EndingDate, averagedealerbalance, OverdraftFund) %>% 
    filter(!is.na(OverdraftFund))

graph=
    ggplot(dayoverdraft4plot1, aes(x = averagedealerbalance, y = OverdraftFund ) ) +
    geom_point(size=4, color = 'blue') +
    geom_smooth(method='lm', formula= y~x, se = FALSE, color = 'black', size = 2) +
    geom_point(data = subset(dayoverdraft4plot1, EndingDate == ymd( '2019-09-25', tz = 'EST' ) ), aes(x = averagedealerbalance, y = OverdraftFund),  size=4, color = 'red') +
    geom_point(data = subset(dayoverdraft4plot1, EndingDate == ymd( '2020-03-11', tz = 'EST' ) ), aes(x = averagedealerbalance, y = OverdraftFund),  size=4, color = 'black', alpha = 1) +
    # coord_cartesian(xlim = c(395,830), ylim = c(0,20), expand = FALSE ) +
    labs(
        # x = TeX('Average daily dealer-bank opening balances ($ billions)'),
        x = ('Average daily top 100 bank opening balances ($ billions)'),
        y = 'Peak daylight overdraft ($ billions)') +
    # scale_x_continuous(  expand = c(0,0)    ) +
    # scale_y_continuous(
    #     # the following specification should not be changed
    #     expand = c(0,0),   # Align plots on y axis
    # ) +      
    # labs(labels = c("Marketable treasuries outstanding", "projected")) +
    # scale_color_discrete(value = c("Legend") ) +
    theme(text = element_text(size = 27),       
          axis.text = element_text(colour = "black", size = rel(0.8)),
          axis.text.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit = "in")   ),
          axis.text.y = element_text(margin = margin(t=0, r=0.1, b=0,l=0, unit = "in")   ),
          # 
          panel.background = element_rect(fill = "white"),
          axis.line = element_line(colour = "black", 
                                   size = 0.1, linetype = "solid"),
          axis.ticks = element_line(colour = "black"),
          axis.ticks.length = unit(-0.05,"in"),
          axis.title.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit='in') ),
          # 
          plot.margin = unit(c(0.5,0.5,0.2,0.6), "in"),  # top right bottom left
          panel.grid.major = element_line(colour = "grey80"),
          # panel.grid.minor = element_line(colour = "grey80")    
          # legend.title =element_blank(), 
          # legend.position= c(0.8, 0.95)
    ) 

graph
# ggsave('OverdraftScatter_top100balances202306.pdf',graph,width=16,height=9, units = "in" )
rm(Data2plot_tem, dayoverdraft)









# ---------------------Time series of payment delay------------------------------------------------------

rm( list = setdiff( ls(), c(keepVariables, "") ))
Data2plot = data2work1 %>% mutate(Date = ymd(Date)) 
data_tem = Data2plot %>% select(Date, QuarterEnd) %>% filter(QuarterEnd>0)

graph1=
    ggplot(Data2plot, aes(x = Date)) +
    geom_line(aes(y = tdiff50_recv1), color='blue',linewidth=0.3) +
    scale_x_date(name = element_blank(),
                     # breaks = date_breaks("1 hour"),
                     breaks=date_breaks('6 months'),
                     # minor_breaks = date_breaks('6 months'),
                     labels = date_format("%m/%Y", tz = "EST"),
                     # limits = c(Regression4$Date[2],Regression4$ Date[328])
                     expand = c(0,0)
    ) +
    # the following specification should not be changed
    scale_y_continuous(name = "Medium Time of Receives (minutes)",
                       expand = c(0,0),   # Align plots on y axis
                       # limits = c(0,27),
                       # n.breaks = 6
    ) +
    geom_vline(xintercept = as.numeric(data_tem$Date), linewidth = 0.4, color = 'orange' ) +
    theme(
        text = element_text(size = 25),
        axis.text = element_text(colour = "black", size = rel(0.8)),
        axis.text.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit = "in"), angle=45, hjust = 1  ),
        axis.text.y = element_text(margin = margin(t=0, r=0.1, b=0,l=0, unit = "in")   ),
        #
        panel.background = element_rect(fill = "white"),
        axis.line = element_line(colour = "black",
                                 size = 0.8, linetype = "solid"),
        axis.ticks = element_line(colour = "black"),
        axis.ticks.length = unit(-0.05,"in"),
        axis.title.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit='in') ),
        #
        plot.margin = unit(c(0.5,0.5,0.2,0.6), "in"),  # top right bottom left
        # panel.grid.major = element_line(colour = "grey80"),
        # panel.grid.minor = element_line(colour = "grey80")
        legend.position = c(0.15, 0.95),
        legend.title = element_blank()
    )
# ggplotly()
graph1
ggsave('MediumTimeSeriesVSQuarterEnd.pdf',graph1,width=16,height=9, units = "in" )






# ---------------------Figure 1------------------------------------------------------
rm( list = setdiff( ls(), c(keepVariables, "") ))
Data2plot = data2work %>% mutate(Date = ymd(Date)) 
cat("\014")  

graphics.off()
Data2plot_tem = Data2plot %>% 
    filter(Date >= ymd('2019-06-15', tz = "EST") ) %>% 
    filter(Date < ymd('2020-11-01', tz = "EST") )
graph1=
    ggplot(Data2plot_tem, aes(x = Date)) +
    # 
    geom_line(aes(y = SOFR_IOER,    color= "SOFR_IOER"   ),  linewidth=1.2) +
    geom_line(aes(y = Repo_balances1*1000/2.5-170,     color= "Repo_balances"    ),  linewidth=0.8) +
    # 
    scale_x_date(
        name = element_blank(),
        # breaks = date_breaks("1 hour"),
        labels = date_format("%d%b%Y", tz = "EST"),
        breaks=date_breaks('1 month'),
        # limits = c(Regression4$Date[2],Regression4$ Date[328])
        expand = c(0,0)
    ) +
    # the following specification should not be changed
    scale_y_continuous(name = "SOFR - IOER (basis points)",
                       expand = c(0,0),   # Align plots on y axis
                       limits = c(-50,350),
                       n.breaks = 9,
                       sec.axis = sec_axis( ~.*2.5+170*2.5, name = "Large dealer balances ($ billion)", breaks = seq(300,1300,100)  )    ) +
    scale_color_manual(name = "Legend",
                       values=c( "SOFR_IOER" = "red", "Repo_balances" = "blue"),
                       labels = c( "SOFR_IOER" = "SOFR - IOER", "Repo_balances" = "Balances")
    ) +
    # scale_color_discrete(name = "Legend") +
    # scale_linetype_manual(name = "Legend",
    #                       values=c( "SOFR_IOER" = "solid", "GCF_IOER" = "longdash", "TGCR_IOER" = "dashed"),
    #                       labels = c( "SOFR_IOER" = "SOFR-IOER", "GCF_IOER" = "GCF-IOER", "TGCR_IOER" = "TGCR-IOER" )
    #                       ) +
    theme(
        text = element_text(size = 20),       
        axis.text = element_text(colour = "black", size = rel(0.8)),
        axis.text.x = element_text(margin = margin(t=0, r=0, b=0, l=0, unit = "in"), angle = 45, hjust = 1, size = 15 ),
        axis.text.y = element_text(margin = margin(t=0, r=0.1, b=0,l=0, unit = "in")   ),
        # axis.text.x = element_text(angle = 45),
        axis.title.y.right = element_text(angle = 90),
        # 
        panel.background = element_rect(fill = "white"),
        axis.line = element_line(colour = "black", 
                                 size = 0.8, linetype = "solid"),
        axis.ticks = element_line(colour = "black"),
        axis.ticks.length = unit(-0.05,"in"),
        axis.title.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit='in') ),
        # 
        plot.margin = unit(c(0.5,0.5,0.2,0.6), "in"),  # top right bottom left
        panel.grid.major = element_line(colour = "grey80"),
        panel.grid.minor = element_line(colour = "grey80"),
        legend.position = c(0.90, 0.95),
        legend.key = element_rect(colour = "transparent", fill = "transparent"), 
        legend.text=element_text(size=15),
        legend.margin = margin(t = -0.5, unit='cm'),
        legend.box.background = element_rect(colour = "black", linewidth = 1),
        legend.box.margin = margin(1, 2, 1, 1),
        legend.title = element_blank()
    ) 
graph1
# ggsave('SOFRvsDealerBalance2.pdf',graph1,width=16,height=9, units = "in" )









# ---------------------Time series of payment delay including Quarter End observations------------------------------------------------------
# Figure 6
rm( list = setdiff( ls(), c(keepVariables, "") ))

data2work = data2work1 %>% 
    # filter(QuarterEnd<1) %>%
    mutate(Ave_tdiff50_recv1_3 = zoo::rollmean(tdiff50_recv1, k = 3, align = "right", fill = NA, na.rm = TRUE),
           Ave_tdiff50_recv1_5 = zoo::rollmean(tdiff50_recv1, k = 5, align = "right", fill = NA, na.rm = TRUE),
           Ave_tdiff50_recv1_7 = zoo::rollmean(tdiff50_recv1, k = 7, align = "right", fill = NA, na.rm = TRUE),
           Ave_tdiff50_recv1_9 = zoo::rollmean(tdiff50_recv1, k = 9, align = "right", fill = NA, na.rm = TRUE),
           Ave_tdiff50_recv1_11 = zoo::rollmean(tdiff50_recv1, k = 11, align = "right", fill = NA, na.rm = TRUE),
           Ave_tdiff50_recv1_13 = zoo::rollmean(tdiff50_recv1, k = 13, align = "right", fill = NA, na.rm = TRUE),
           Ave_tdiff50_recv1_15 = zoo::rollmean(tdiff50_recv1, k = 15, align = "right", fill = NA, na.rm = TRUE),
           Ave_tdiff50_recv1_17 = zoo::rollmean(tdiff50_recv1, k = 17, align = "right", fill = NA, na.rm = TRUE),
           Ave_tdiff50_recv1_19 = zoo::rollmean(tdiff50_recv1, k = 19, align = "right", fill = NA, na.rm = TRUE),
           Ave_tdiff50_recv1_21 = zoo::rollmean(tdiff50_recv1, k = 21, align = "right", fill = NA, na.rm = TRUE),
           AveBalancesPreviousWeeklag3 = lag(Ave_tdiff50_recv1_3, 2),
           AveBalancesPreviousWeeklag5 = lag(Ave_tdiff50_recv1_5, 2),
           AveBalancesPreviousWeeklag7 = lag(Ave_tdiff50_recv1_7, 2),
           AveBalancesPreviousWeeklag9 = lag(Ave_tdiff50_recv1_9, 2),
           AveBalancesPreviousWeeklag11 = lag(Ave_tdiff50_recv1_11, 2),
           AveBalancesPreviousWeeklag13 = lag(Ave_tdiff50_recv1_13, 2),
           AveBalancesPreviousWeeklag15 = lag(Ave_tdiff50_recv1_15, 2),
           AveBalancesPreviousWeeklag17 = lag(Ave_tdiff50_recv1_17, 2),
           AveBalancesPreviousWeeklag19 = lag(Ave_tdiff50_recv1_19, 2),
           AveBalancesPreviousWeeklag21 = lag(Ave_tdiff50_recv1_21, 2),
           lag2_Repo_balances1 = lag(Repo_balances1, 2)
    )

# mean(data2work$tdiff50_recv1, na.rm = TRUE)
# mean(data2work1$tdiff50_recv1, na.rm = TRUE)
# mean(data2work$Ave_tdiff50_recv1_11, na.rm = TRUE)
# data2work2 = data2work1 %>% 
#     filter(QuarterEnd>=1) %>% select(Date, tdiff50_recv1)
# mean(data2work2$tdiff50_recv1, na.rm = TRUE)

Data2plot = data2work %>% mutate(Date = ymd(Date)) %>% filter(!is.na(Ave_tdiff50_recv1_9))

graph1=
    ggplot(Data2plot, aes(x = Date)) +
    geom_rect(aes( xmin= as.Date('2020-03-01') , xmax = as.Date('2020-04-30') , ymin=-Inf, ymax=Inf ), alpha=0.01, fill = 'pink') +
    geom_line(aes(y = Ave_tdiff50_recv1_9), color='blue',linewidth=1) +
    geom_vline(xintercept = as.numeric( ymd('2019-09-18') ), linewidth = 0.8, color = 'orange' ) +
    # the following specification should not be changed
    scale_y_continuous(name = "Lagged average payment delay (minutes)",
                       expand = c(0,0),   # Align plots on y axis
                       limits = c( min(Data2plot$Ave_tdiff50_recv1_11)-10 ,  max(Data2plot$Ave_tdiff50_recv1_11)+10),
                       # n.breaks = 6
    ) +
    scale_x_date(name = element_blank(),
                 # breaks = date_breaks("1 hour"),
                 breaks=date_breaks('6 months'),
                 # minor_breaks = date_breaks('6 months'),
                 labels = date_format("%m/%Y"),
                 # limits = c(Regression4$Date[2],Regression4$ Date[328])
                 expand = c(0,0)
    ) +
    theme(
        text = element_text(size = 25),
        axis.text = element_text(colour = "black", size = rel(0.8)),
        axis.text.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit = "in"), angle=45, hjust = 1  ),
        axis.text.y = element_text(margin = margin(t=0, r=0.1, b=0,l=0, unit = "in")   ),
        #
        panel.background = element_rect(fill = "white"),
        axis.line = element_line(colour = "black",
                                 linewidth = 0.8, linetype = "solid"),
        axis.ticks = element_line(colour = "black"),
        axis.ticks.length = unit(-0.05,"in"),
        axis.title.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit='in') ),
        #
        plot.margin = unit(c(0.5,0.5,0.2,0.6), "in"),  # top right bottom left
        # panel.grid.major = element_line(colour = "grey80"),
        # panel.grid.minor = element_line(colour = "grey80")
        legend.position = c(0.15, 0.95),
        legend.title = element_blank()
    )
# ggplotly()
graph1
# ggsave('MediumTimeSeriesAverage9BusinessDaysIncludingQuarterEnd.pdf',graph1,width=16,height=9, units = "in" )

# ggsave('DelayPlot.pdf',graph1,width=16,height=9, units = "in" )

min(Data2plot$Ave_tdiff50_recv1_9, na.rm = TRUE)

# tem = Data2plot %>% select(Date, Ave_tdiff50_recv1_9)



# ---------------------Time series of payment delay including Quarter End observations------------------------------------------------------

rm( list = setdiff( ls(), c(keepVariables, "") ))


data2worktem = data2work1 %>% 
    mutate(
        tdiff50_recv1 = if_else(QuarterEnd>=1, NA, tdiff50_recv1)
    ) %>% select(Date, tdiff50_recv1) %>%
    mutate(Ave_tdiff50_recv1_3 = zoo::rollmean(tdiff50_recv1, k = 3, align = "right", fill = NA, na.rm = TRUE),
           Ave_tdiff50_recv1_5 = zoo::rollmean(tdiff50_recv1, k = 5, align = "right", fill = NA, na.rm = TRUE),
           Ave_tdiff50_recv1_7 = zoo::rollmean(tdiff50_recv1, k = 7, align = "right", fill = NA, na.rm = TRUE),
           Ave_tdiff50_recv1_9 = zoo::rollmean(tdiff50_recv1, k = 9, align = "right", fill = NA, na.rm = TRUE),
           Ave_tdiff50_recv1_11 = zoo::rollmean(tdiff50_recv1, k = 11, align = "right", fill = NA, na.rm = TRUE),
           Ave_tdiff50_recv1_13 = zoo::rollmean(tdiff50_recv1, k = 13, align = "right", fill = NA, na.rm = TRUE),
           Ave_tdiff50_recv1_15 = zoo::rollmean(tdiff50_recv1, k = 15, align = "right", fill = NA, na.rm = TRUE),
           Ave_tdiff50_recv1_17 = zoo::rollmean(tdiff50_recv1, k = 17, align = "right", fill = NA, na.rm = TRUE),
           Ave_tdiff50_recv1_19 = zoo::rollmean(tdiff50_recv1, k = 19, align = "right", fill = NA, na.rm = TRUE),
           Ave_tdiff50_recv1_21 = zoo::rollmean(tdiff50_recv1, k = 21, align = "right", fill = NA, na.rm = TRUE),
           AveBalancesPreviousWeeklag3 = lag(Ave_tdiff50_recv1_3, 2),
           AveBalancesPreviousWeeklag5 = lag(Ave_tdiff50_recv1_5, 2),
           AveBalancesPreviousWeeklag7 = lag(Ave_tdiff50_recv1_7, 2),
           AveBalancesPreviousWeeklag9 = lag(Ave_tdiff50_recv1_9, 2),
           AveBalancesPreviousWeeklag11 = lag(Ave_tdiff50_recv1_11, 2),
           AveBalancesPreviousWeeklag13 = lag(Ave_tdiff50_recv1_13, 2),
           AveBalancesPreviousWeeklag15 = lag(Ave_tdiff50_recv1_15, 2),
           AveBalancesPreviousWeeklag17 = lag(Ave_tdiff50_recv1_17, 2),
           AveBalancesPreviousWeeklag19 = lag(Ave_tdiff50_recv1_19, 2),
           AveBalancesPreviousWeeklag21 = lag(Ave_tdiff50_recv1_21, 2)) %>% 
    select(Date, 
           Ave_tdiff50_recv1_3 ,
           Ave_tdiff50_recv1_5 ,
           Ave_tdiff50_recv1_7 ,
           Ave_tdiff50_recv1_9 ,
           Ave_tdiff50_recv1_11,
           Ave_tdiff50_recv1_13,
           Ave_tdiff50_recv1_15,
           Ave_tdiff50_recv1_17,
           Ave_tdiff50_recv1_19,
           Ave_tdiff50_recv1_21,
           AveBalancesPreviousWeeklag3 ,
           AveBalancesPreviousWeeklag5 ,
           AveBalancesPreviousWeeklag7 ,
           AveBalancesPreviousWeeklag9 ,
           AveBalancesPreviousWeeklag11,
           AveBalancesPreviousWeeklag13,
           AveBalancesPreviousWeeklag15,
           AveBalancesPreviousWeeklag17,
           AveBalancesPreviousWeeklag19,
           AveBalancesPreviousWeeklag21
    )

data2work = left_join(data2work1, data2worktem, by = 'Date')

# mean(data2work$tdiff50_recv1, na.rm = TRUE)
# mean(data2work1$tdiff50_recv1, na.rm = TRUE)
# mean(data2work$Ave_tdiff50_recv1_11, na.rm = TRUE)
# data2work2 = data2work1 %>% 
#     filter(QuarterEnd>=1) %>% select(Date, tdiff50_recv1)
# mean(data2work2$tdiff50_recv1, na.rm = TRUE)

Data2plot = data2work %>% mutate(Date = ymd(Date)) %>% filter(!is.na(Ave_tdiff50_recv1_9))

graph1=
    ggplot(Data2plot, aes(x = Date)) +
    geom_rect(aes( xmin= as.Date('2020-03-01') , xmax = as.Date('2020-04-30') , ymin=-Inf, ymax=Inf ), alpha=0.01, fill = 'pink') +
    geom_line(aes(y = Ave_tdiff50_recv1_9), color='blue',linewidth=1) +
    geom_vline(xintercept = as.numeric( ymd('2019-09-18') ), linewidth = 0.8, color = 'orange' ) +
    # the following specification should not be changed
    scale_y_continuous(name = "Lagged average payment delay (minutes)",
                       expand = c(0,0),   # Align plots on y axis
                       limits = c( min(Data2plot$Ave_tdiff50_recv1_11)-10 ,  max(Data2plot$Ave_tdiff50_recv1_11)+10),
                       # n.breaks = 6
    ) +
    scale_x_date(name = element_blank(),
                 # breaks = date_breaks("1 hour"),
                 breaks=date_breaks('6 months'),
                 # minor_breaks = date_breaks('6 months'),
                 labels = date_format("%m/%Y"),
                 # limits = c(Regression4$Date[2],Regression4$ Date[328])
                 expand = c(0,0)
    ) +
    theme(
        text = element_text(size = 25),
        axis.text = element_text(colour = "black", size = rel(0.8)),
        axis.text.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit = "in"), angle=45, hjust = 1  ),
        axis.text.y = element_text(margin = margin(t=0, r=0.1, b=0,l=0, unit = "in")   ),
        #
        panel.background = element_rect(fill = "white"),
        axis.line = element_line(colour = "black",
                                 linewidth = 0.8, linetype = "solid"),
        axis.ticks = element_line(colour = "black"),
        axis.ticks.length = unit(-0.05,"in"),
        axis.title.x = element_text(margin = margin(t=0.1, r=0, b=0, l=0, unit='in') ),
        #
        plot.margin = unit(c(0.5,0.5,0.2,0.6), "in"),  # top right bottom left
        # panel.grid.major = element_line(colour = "grey80"),
        # panel.grid.minor = element_line(colour = "grey80")
        legend.position = c(0.15, 0.95),
        legend.title = element_blank()
    )
# ggplotly()
graph1
# ggsave('MediumTimeSeriesAverage9BusinessDaysExcludingQuarterEnd.pdf',graph1,width=16,height=9, units = "in" )



